Segala puji dan syukur penulis panjatkan ke hadirat Tuhan Yang Maha Esa atas limpahan rahmat, karunia, dan kemudahan-Nya sehingga eBook ini yang berjudul “Analisis Data Kategori” dapat diselesaikan dengan baik.
Penulisan ini disusun sebagai bagian dari pemenuhan tugas mata kuliah Analisis Data Kategori yang diampu oleh Dr. I Gede Nyoman Mindra Jaya pada semester genap tahun akademik 2024/2025, di Program Studi Statistika, FMIPA Universitas Padjadjaran.
Sebagai mahasiswa semester empat, penulis menyadari bahwa pemahaman terhadap analisis data kategori sangat penting, terutama karena metode-metode seperti tabel kontingensi, uji chi-square, regresi logistik, hingga model log-linear dan Generalized Linear Model (GLM) merupakan materi inti yang sering digunakan dalam berbagai penelitian sosial, kesehatan, bisnis, maupun bidang statistik lainnya.
Oleh karena itu, penyusunan eBook ini tidak hanya ditujukan untuk memenuhi kewajiban akademik, tetapi juga sebagai sarana pembelajaran mandiri yang diharapkan dapat membantu rekan-rekan mahasiswa dalam memahami materi secara lebih utuh dan aplikatif.
Dengan mengacu pada struktur pembelajaran yang disampaikan selama perkuliahan serta referensi yang diberikan oleh dosen pengampu, penulis berupaya merangkum isi buku ini secara sistematis dan ringkas, namun tetap memuat contoh-contoh praktis menggunakan perangkat lunak R agar lebih relevan dan mudah dipahami.
Harapannya, eBook ini dapat menjadi bahan belajar tambahan, terutama untuk teman-teman yang masih merasa kesulitan dalam menginterpretasikan hasil-hasil analisis data kategori.
Ucapan terima kasih yang tulus penulis sampaikan kepada Dr. I Gede Nyoman Mindra Jaya, selaku dosen pengampu mata kuliah, atas ilmu, penjelasan yang runut, serta bimbingan dan motivasi yang telah diberikan selama proses pembelajaran berlangsung.
Penulis juga mengapresiasi dukungan dan diskusi yang berharga dari rekan-rekan sekelas yang ikut berkontribusi dalam memperluas pemahaman terhadap materi.
Penulis menyadari bahwa eBook ini masih memiliki berbagai kekurangan, baik dari segi penyusunan maupun kedalaman materi. Oleh karena itu, masukan, kritik, dan saran yang membangun sangat penulis harapkan untuk penyempurnaan di masa mendatang.
Akhir kata, semoga eBook ini dapat memberikan manfaat nyata dan menjadi referensi yang berguna bagi siapa pun yang ingin mempelajari lebih lanjut tentang analisis data kategori.
Dalam analisis statistik, data kategori merupakan tipe data yang mengelompokkan objek ke dalam kelas-kelas tertentu berdasarkan karakteristik atau atribut yang dimiliki. Tidak seperti data numerik yang bersifat kontinu, data kategori tidak memiliki nilai urutan yang tetap atau jarak antar kategori yang terukur. Beberapa contoh data kategori yang umum digunakan antara lain adalah jenis kelamin (pria/wanita), status pekerjaan (bekerja/tidak bekerja), tingkat kepuasan pelanggan (puas/tidak puas), dan pilihan politik.
Data kategori sering menjadi fokus dalam berbagai studi karena kemampuannya merepresentasikan realitas sosial dan perilaku manusia dalam bentuk klasifikasi. Untuk menganalisis data seperti ini, pendekatan yang digunakan biasanya berbeda dari analisis data numerik. Teknik-teknik seperti tabel kontingensi, uji independensi (uji chi-square), serta model probabilistik seperti regresi logistik, digunakan untuk mengidentifikasi keterkaitan antar kategori.
Dengan kemajuan teknologi informasi, analisis data kategori juga mulai diintegrasikan dalam sistem kecerdasan buatan dan algoritma machine learning. Dalam banyak kasus, data ini harus diubah ke dalam format numerik terlebih dahulu, menggunakan teknik seperti label encoding atau one-hot encoding, agar dapat diproses lebih lanjut. Analisis data kategori tidak hanya memberikan pemahaman deskriptif, tetapi juga berperan dalam membangun sistem prediktif dan pengambilan keputusan berbasis data.
Analisis terhadap data kategori bertujuan untuk menggali wawasan yang tersembunyi dari data yang bersifat klasifikatif. Tujuan ini mencakup beberapa aspek penting berikut:
Melalui pengelompokan dan tabulasi data kategori, peneliti dapat mendeteksi pola-pola umum yang muncul dalam kelompok tertentu. Misalnya, dalam bidang pemasaran, preferensi konsumen terhadap suatu produk dapat dianalisis berdasarkan usia, jenis kelamin, atau wilayah domisili.
Struktur Matematis: Pendekatan ini melibatkan penggunaan frekuensi relatif, persentase, dan representasi visual seperti diagram batang atau pie chart. Distribusi probabilitas kategori juga sering dihitung untuk melihat proporsi antar pilihan.
Penelitian yang melibatkan dua atau lebih variabel kategori sering mencari tahu apakah ada hubungan signifikan antar variabel tersebut. Misalnya, apakah terdapat kaitan antara latar belakang pendidikan dengan pilihan karier, atau antara gaya hidup dengan kesehatan.
Teknik Analisis: Pengujian statistik seperti chi-square digunakan untuk mengevaluasi apakah hubungan yang tampak dalam data terjadi secara kebetulan atau memiliki signifikansi statistik. Untuk analisis lebih lanjut, digunakan juga regresi logistik apabila salah satu variabel merupakan variabel dependen kategori.
Temuan dari analisis data kategori bisa dimanfaatkan untuk menyusun strategi atau kebijakan yang lebih efektif. Dalam konteks manajemen publik, informasi ini bisa membantu merancang program berdasarkan klasifikasi wilayah atau kelompok masyarakat tertentu.
Dengan pendekatan berbasis bukti ini, intervensi dapat diarahkan secara lebih terarah dan efisien, mengingat strategi disesuaikan dengan karakteristik kategori yang dominan dalam populasi.
Dalam dunia modern yang didukung oleh kecerdasan buatan, data kategori sering kali menjadi bagian dari input model prediksi. Misalnya, prediksi terhadap keputusan pembelian pelanggan berdasarkan kategori preferensi, jenis kelamin, dan histori transaksi.
Agar dapat diproses oleh algoritma machine learning, data kategori diubah menjadi format numerik menggunakan metode seperti one-hot encoding. Setelah diolah, model prediktif seperti regresi logistik atau decision tree dapat digunakan untuk memperkirakan kemungkinan suatu kejadian berdasarkan kombinasi dari beberapa variabel kategori.
Analisis ini digunakan secara luas di berbagai bidang:
Analisis data kategori adalah analisis statistik terhadap variabel yang berupa kategori, bukan nilai numerik kontinu.
Kasus: Perusahaan ingin mengevaluasi tingkat
kepuasan melalui komentar sosial media.
Kategori: Sangat Puas, Puas, Netral, Tidak Puas, Sangat
Tidak Puas
Metode: Uji Chi-Square untuk melihat hubungan antara
sentimen dan faktor usia.
Kasus: Bank ingin memprediksi loyalitas berdasarkan
usia dan jenis akun.
Metode: Regresi logistik biner.
Kasus: Rumah sakit mengkaji hubungan antara merokok,
tekanan darah, jenis kelamin, dan penyakit jantung.
Metode: Model loglinear.
Kasus: Mengelompokkan mahasiswa berdasarkan nilai,
kehadiran, dan kebiasaan belajar.
Metode: K-means / k-modes clustering.
Analisis data kategori memiliki beragam pendekatan tergantung pada struktur data dan tujuan analisisnya. Metode-metode ini tidak hanya berfungsi untuk menjelaskan pola dan hubungan, tetapi juga dapat dimanfaatkan untuk membangun model prediksi dan visualisasi eksploratif. Dalam bagian ini akan dibahas empat metode utama yang sering digunakan, dimulai dari teknik dasar seperti tabel kontingensi hingga pendekatan machine learning seperti random forest.
Tabel kontingensi adalah alat penting dalam analisis data kategori karena menyajikan frekuensi observasi berdasarkan kombinasi dua atau lebih variabel kategori. Tabel ini menjadi dasar bagi berbagai pengujian statistik, terutama uji Chi-Square.
Sebagai contoh sederhana:
| Sakit | Tidak Sakit | |
|---|---|---|
| Merokok | 30 | 10 |
| Tidak Merokok | 15 | 45 |
Tabel ini membantu kita mengamati distribusi frekuensi dan memeriksa apakah distribusi tersebut terjadi secara acak atau menunjukkan pola yang signifikan.
Uji Chi-Square digunakan untuk mengevaluasi apakah ada hubungan yang signifikan antara dua variabel kategori. Misalnya, untuk mengetahui apakah status merokok berkaitan dengan kondisi kesehatan.
# Contoh uji chi-square di R
data <- matrix(c(30,10,15,45), nrow = 2, byrow = TRUE)
colnames(data) <- c("Sakit", "Tidak Sakit")
rownames(data) <- c("Merokok", "Tidak Merokok")
chisq.test(data)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: data
## X-squared = 22.264, df = 1, p-value = 2.376e-06
Interpretasi:
Regresi logistik merupakan teknik yang digunakan ketika variabel dependen bersifat biner, seperti ya/tidak atau sukses/gagal. Pendekatan ini dapat menangani prediktor kategori maupun numerik.
set.seed(123)
data <- data.frame(
loyal = sample(0:1, 100, replace = TRUE),
usia = sample(18:60, 100, replace = TRUE),
frek_transaksi = rpois(100, lambda = 5)
)
# Regresi logistik
model <- glm(loyal ~ usia + frek_transaksi, data = data, family = "binomial")
summary(model)
Interpretasi Koefisien:
exp(koefisien).Regresi logistik menjadi sangat berguna dalam analisis prediktif berbasis kategori, seperti menentukan kemungkinan pelanggan melakukan pembelian atau tidak.
Correspondence Analysis adalah teknik eksplorasi yang digunakan untuk menganalisis dan memvisualisasikan hubungan antara kategori dalam tabel kontingensi yang lebih besar. CA memiliki peran serupa dengan PCA (Principal Component Analysis), namun dirancang untuk data kategori.
library(FactoMineR)
library(factoextra)
data <- matrix(c(20,30,50,40,60,30,25,45,35), nrow=3, byrow=TRUE)
rownames(data) <- c("Remaja", "Dewasa", "Lansia")
colnames(data) <- c("Makanan A", "Makanan B", "Makanan C")
res.ca <- CA(data, graph = FALSE)
fviz_ca_biplot(res.ca)
Dengan CA, kita dapat mengamati asosiasi visual antara preferensi makanan dan kelompok usia dalam ruang dua dimensi.
Metode pohon keputusan (decision tree) dan hutan acak (random forest) merupakan teknik klasifikasi berbasis machine learning yang populer untuk data kategori. Decision tree memberikan aturan klasifikasi yang mudah dipahami, sedangkan random forest memperbaiki akurasi melalui pendekatan ensemble.
library(rpart)
set.seed(123)
data <- data.frame(
loyal = factor(sample(c("ya", "tidak"), 100, replace = TRUE)),
usia = sample(18:60, 100, replace = TRUE),
jenis_akun = factor(sample(c("Premium", "Reguler"), 100, replace = TRUE))
)
tree <- rpart(loyal ~ usia + jenis_akun, data = data, method = "class")
library(rpart.plot)
rpart.plot(tree)
library(randomForest)
rf <- randomForest(loyal ~ ., data = data, ntree = 100)
print(rf)
| Metode | Kelebihan | Kekurangan |
|---|---|---|
| Decision Tree | Interpretasi mudah | Cenderung overfitting |
| Random Forest | Akurasi tinggi, stabil | Interpretasi rumit |
Secara keseluruhan, keempat metode ini menawarkan pendekatan yang saling melengkapi dalam menganalisis data kategori. Tabel kontingensi dan uji chi-square unggul untuk asosiasi dasar, regresi logistik kuat dalam prediksi biner, correspondence analysis menyediakan eksplorasi visual, dan metode machine learning seperti decision tree dan random forest menawarkan solusi klasifikasi berbasis kategori yang adaptif. Kombinasi dari metode-metode ini dapat memberikan pemahaman yang komprehensif terhadap struktur data kategori yang kompleks.
Distribusi probabilitas berperan penting dalam pemodelan kejadian acak yang bersifat kategori atau diskrit. Dalam bagian ini, akan dibahas empat jenis distribusi probabilitas yang sering digunakan dalam analisis data kategori: Bernoulli, Binomial, Multinomial, dan Poisson. Penjelasan mencakup bentuk fungsi probabilitas, interpretasi notasi matematis, serta contoh implementasi di R.
Distribusi Bernoulli digunakan untuk eksperimen dengan dua kemungkinan hasil: sukses (1) dan gagal (0).
Fungsi probabilitas:
\(P(X = x) = p^x (1 - p)^{1 - x}, \quad x \in
\{0, 1\}\)
Contoh dalam R:
set.seed(123)
rbinom(n = 10, size = 1, prob = 0.5)
## [1] 0 1 0 1 1 0 1 1 1 0
Contoh kasus: Menentukan apakah seorang siswa lulus (1) atau tidak lulus (0) ujian.
Distribusi Binomial adalah perpanjangan dari Bernoulli untuk \(n\) percobaan independen.
Fungsi probabilitas:
\(P(X = k) = \binom{n}{k} p^k
(1-p)^{n-k}\)
Contoh dalam R:
set.seed(123)
rbinom(n = 10, size = 5, prob = 0.5)
## [1] 2 3 2 4 4 1 3 4 3 2
Contoh kasus: Berapa banyak dari 5 mahasiswa yang lulus ujian jika tiap mahasiswa punya peluang 50%.
Distribusi Multinomial merupakan generalisasi distribusi Binomial untuk lebih dari dua kategori.
Fungsi probabilitas:
\(P(X_1 = x_1, ..., X_k = x_k) =
\frac{n!}{x_1!x_2!...x_k!} p_1^{x_1} p_2^{x_2} ...
p_k^{x_k}\)
Contoh dalam R:
set.seed(123)
rmultinom(n = 1, size = 10, prob = c(0.3, 0.5, 0.2))
## [,1]
## [1,] 2
## [2,] 5
## [3,] 3
Contoh kasus: Memilih 10 orang dan mencatat pilihan makanan favorit mereka dari 3 jenis.
Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam suatu interval waktu atau ruang tertentu, terutama jika kejadian tersebut jarang terjadi.
Fungsi probabilitas:
\(P(X = k) = \frac{e^{-\lambda}
\lambda^k}{k!}\)
Contoh dalam R:
set.seed(123)
rpois(10, lambda = 3)
## [1] 2 4 2 5 6 0 3 5 3 3
Contoh kasus: Jumlah telepon masuk ke pusat layanan pelanggan dalam satu jam.
Distribusi-distribusi di atas merupakan dasar yang sangat penting dalam menganalisis dan memodelkan fenomena kategorik. Pemahaman terhadap bentuk dan notasinya akan membantu dalam memilih model yang sesuai dan melakukan interpretasi hasil analisis secara tepat.
Dalam analisis data kategori, pemilihan desain sampling merupakan tahapan penting yang menentukan kualitas hasil penelitian, baik dari segi validitas maupun reliabilitasnya. Desain sampling yang dipilih harus selaras dengan tujuan riset serta jenis data yang hendak dikumpulkan. Secara garis besar, desain sampling dalam konteks ini terbagi menjadi dua pendekatan utama, yakni prospective sampling dan retrospective sampling.
Kedua pendekatan tersebut memiliki karakteristik tersendiri yang mempengaruhi metode sampling yang digunakan. Dalam praktiknya, desain sampling ini umum dijumpai dalam berbagai jenis penelitian, baik yang bersifat eksperimental maupun observasional, seperti eksperimen, studi kohort, dan studi kasus-kontrol. Pemahaman terhadap kerangka desain sampling ini menjadi kunci dalam merancang penelitian yang efektif, efisien, dan mampu memberikan bukti yang meyakinkan.
Prospective sampling merupakan pendekatan di mana subjek penelitian dipilih dan diikuti secara berkelanjutan sepanjang periode waktu tertentu. Tujuan utamanya adalah untuk mencatat atau mengamati perkembangan variabel yang diteliti secara langsung. Dengan pendekatan ini, peneliti memiliki kesempatan untuk mengendalikan variabel bebas sebelum proses pengukuran dilakukan. Oleh karena itu, metode ini sangat cocok untuk studi yang ingin menguji hubungan kausal, terutama dalam penelitian eksperimental.
Studi eksperimental adalah jenis penelitian yang dirancang untuk menguji pengaruh suatu perlakuan (treatment) dengan mengontrol variabel lain. Dalam pendekatan ini, subjek dialokasikan secara acak ke dalam kelompok perlakuan dan kontrol. Teknik sampling yang umum digunakan meliputi:
Langkah pelaksanaan:
Konteks penggunaan:
Digunakan ketika penelitian bertujuan untuk menemukan hubungan sebab-akibat secara langsung dengan intervensi yang dapat dikontrol.
Studi kohort merupakan bentuk penelitian observasional prospektif yang mengikuti sekelompok orang berdasarkan status paparan mereka terhadap suatu faktor risiko. Subjek dengan dan tanpa paparan diikuti sepanjang waktu untuk mencatat apakah mereka mengembangkan outcome tertentu.
Metode sampling yang digunakan:
Langkah pelaksanaan:
Konteks penggunaan:
Cocok ketika peneliti ingin mengamati efek jangka panjang dari paparan terhadap outcome yang tidak langsung terlihat.
Pendekatan retrospective sampling menggunakan data yang dikumpulkan dari peristiwa atau kejadian yang sudah terjadi. Peneliti mulai dari outcome terlebih dahulu, kemudian menelusuri variabel paparan masa lalu. Karena efisiensi waktu dan biaya, pendekatan ini sering digunakan dalam studi yang berkaitan dengan kondisi langka atau outcome yang memerlukan waktu lama untuk muncul.
Studi ini dimulai dengan mengidentifikasi individu yang memiliki kondisi atau outcome (kasus) dan membandingkannya dengan individu yang tidak memiliki kondisi tersebut (kontrol). Peneliti kemudian menelusuri informasi masa lalu untuk menentukan apakah paparan terhadap faktor risiko berbeda antara kedua kelompok.
Metode sampling yang digunakan:
Langkah pelaksanaan:
Konteks penggunaan:
Ideal untuk penyakit atau kondisi langka, serta situasi ketika outcome sudah terjadi sebelum penelitian dimulai.
Pendekatan ini menyerupai studi kohort, tetapi dilakukan dengan menggunakan data historis. Peneliti mengelompokkan subjek berdasarkan paparan di masa lalu, kemudian melihat outcome yang sudah terjadi.
Metode sampling yang digunakan:
Konteks penggunaan:
Digunakan ketika informasi historis tersedia dan ingin dianalisis untuk menggambarkan pola hubungan antara paparan dan outcome di masa lalu.
| Jenis Studi | Pendekatan | Metode Sampling | Keuntungan | Kelemahan |
|---|---|---|---|---|
| Eksperimen | Prospektif | SRS, Stratified, Cluster | Kontrol tinggi, dapat menganalisis kausal | Biaya tinggi, isu etika |
| Studi Kohort | Prospektif | Census, Systematic, Matched | Dapat melihat perkembangan jangka panjang | Waktu lama, risiko kehilangan data |
| Studi Kasus-Kontrol | Retrospektif | Purposive, Snowball, Incidence Density | Cepat, efisien untuk penyakit langka | Rentan bias recall |
| Kohort Retrospektif | Retrospektif | Convenience, Quota, Case-Based | Gunakan data historis, biaya rendah | Kualitas data tergantung rekaman |
| Aspek | Eksperimen | Studi Kohort |
|---|---|---|
| Desain | Intervensi | Observasional |
| Randomisasi | Ya | Tidak |
| Biaya & Waktu | Tinggi, lebih cepat | Lebih murah, lebih lama |
| Tujuan | Menilai efek perlakuan | Mengamati efek paparan alami |
| Kontrol Variabel | Tinggi | Terbatas |
| Aspek | Studi Kohort | Studi Kasus-Kontrol |
|---|---|---|
| Arah Pengamatan | Maju (prospektif) atau mundur (retrospektif) | Mundur dari outcome ke paparan |
| Waktu dan Biaya | Lama dan mahal | Cepat dan murah |
| Outcome yang dikaji | Umum | Langka |
| Ukuran yang dihitung | Risk Ratio, Incidence Rate | Odds Ratio |
| Bias | Rendah | Rentan bias recall |
Dengan memahami berbagai pendekatan dan metode sampling yang tersedia, peneliti dapat lebih cermat dalam memilih strategi pengambilan sampel yang sesuai dengan jenis data kategori dan desain penelitian yang digunakan. Pilihan yang tepat akan membantu meningkatkan keakuratan hasil analisis dan validitas kesimpulan penelitian.
Tabel kontingensi merupakan alat utama dalam statistik kategori yang digunakan untuk menyajikan data dalam bentuk matriks frekuensi. Fungsinya adalah untuk menggambarkan hubungan antara dua atau lebih variabel kategori dan menjadi dasar dari berbagai teknik analisis inferensial. Pemahaman mengenai tipe tabel kontingensi, distribusi frekuensi, dan proporsi sangat penting sebelum menerapkan uji statistik seperti chi-square atau model loglinear.
Jenis tabel kontingensi ditentukan berdasarkan jumlah variabel dan kategori yang terlibat. Setiap bentuk tabel menyajikan kompleksitas dan kegunaan yang berbeda sesuai tujuan analisis.
Tabel jenis ini melibatkan dua variabel kategori, masing-masing dengan dua tingkat. Sering digunakan untuk mengevaluasi hubungan sederhana seperti efek suatu kondisi terhadap hasil tertentu.
Contoh:
| Penyakit Paru-paru: Ya | Penyakit Paru-paru: Tidak | |
|---|---|---|
| Merokok: Ya | 30 | 10 |
| Merokok: Tidak | 15 | 45 |
Interpretasi awal dari tabel ini dapat dilakukan dengan membandingkan proporsi di setiap baris dan kolom untuk melihat adanya ketimpangan distribusi.
Tabel RxC memuat dua variabel dengan lebih dari dua kategori. Biasanya digunakan untuk analisis stratifikasi atau efek multikategori terhadap suatu variabel kualitatif.
Contoh:
| Pendidikan | Pengangguran | Pekerja | Wiraswasta |
|---|---|---|---|
| SD | 10 | 20 | 5 |
| SMP | 12 | 25 | 8 |
| SMA | 15 | 30 | 12 |
| S1 | 5 | 40 | 20 |
| S2 | 2 | 18 | 10 |
Dari tabel ini, dapat dilakukan analisis untuk mengidentifikasi apakah tingkat pendidikan berkorelasi dengan jenis pekerjaan.
Jenis ini digunakan ketika terdapat tiga atau lebih variabel kategori. Biasanya ditampilkan dalam bentuk daftar berlapis atau tabel tiga dimensi dan digunakan untuk eksplorasi interaksi kompleks.
Contoh kasus: Hubungan antara Jenis Kelamin, Status Pekerjaan, dan Tingkat Pendidikan. Model loglinear atau analisis interaksi banyak variabel sangat cocok untuk jenis tabel ini.
Setelah menyusun tabel kontingensi, analisis dilanjutkan dengan melihat berbagai jenis distribusi untuk memahami struktur data.
Frekuensi absolut menunjukkan jumlah pengamatan aktual dalam tiap sel. Ini adalah nilai mentah dari data.
Contoh: Sel [Merokok = Ya, Penyakit = Ya] berisi 30 → artinya ada 30 orang dengan karakteristik tersebut.
Frekuensi relatif menyatakan proporsi suatu sel terhadap total observasi.
\(\text{Proporsi} = \frac{\text{Frekuensi Sel}}{\text{Total Observasi}}\)
Contoh dalam R:
tab <- matrix(c(30,10,15,45), nrow = 2, byrow = TRUE)
proporsi <- prop.table(tab)
proporsi
## [,1] [,2]
## [1,] 0.30 0.10
## [2,] 0.15 0.45
Distribusi marginal adalah total frekuensi untuk setiap baris atau kolom, berguna untuk melihat distribusi total dari satu variabel.
Contoh dalam R:
marginal_baris <- margin.table(tab, 1)
marginal_kolom <- margin.table(tab, 2)
marginal_baris
## [1] 40 60
marginal_kolom
## [1] 45 55
Distribusi kondisional menghitung probabilitas dengan kondisi tetap pada baris atau kolom tertentu.
Contoh: Probabilitas terkena penyakit paru-paru jika diketahui status merokok.
cond_row <- prop.table(tab, 1) # per baris
cond_col <- prop.table(tab, 2) # per kolom
cond_row
## [,1] [,2]
## [1,] 0.75 0.25
## [2,] 0.25 0.75
cond_col
## [,1] [,2]
## [1,] 0.6666667 0.1818182
## [2,] 0.3333333 0.8181818
Interpretasi dari proporsi baris dapat mengindikasikan apakah merokok berpengaruh terhadap probabilitas memiliki penyakit paru-paru.
Visualisasi membantu memperjelas pola hubungan antar variabel dalam tabel kontingensi. Salah satu alat visual yang efektif adalah mosaic plot.
mosaicplot(tab, main = "Mosaic Plot Tabel Kontingensi", color = TRUE)
Mosaic plot memperlihatkan proporsi setiap kategori dalam bentuk area, sehingga perbandingan antar sel atau antar kelompok lebih intuitif.
Tabel kontingensi adalah pondasi penting dalam eksplorasi dan analisis data kategori. Memahami bentuk-bentuk tabel dan bagaimana membacanya memungkinkan peneliti untuk mengambil langkah analitik yang tepat dalam tahap selanjutnya seperti uji chi-square, model loglinear, atau regresi logistik. Pengetahuan ini tidak hanya membantu dalam interpretasi data, tetapi juga dalam menyusun argumen statistik berbasis bukti.
Tabel kontingensi 2x2 merupakan bentuk paling dasar dari representasi hubungan dua variabel kategori yang masing-masing memiliki dua kategori. Meskipun sederhana, tabel ini menjadi landasan penting untuk banyak metode statistik inferensial, seperti uji chi-square dan analisis regresi logistik biner. Tabel 2x2 tidak hanya menyajikan frekuensi absolut, tetapi juga memungkinkan perhitungan distribusi peluang berupa probabilitas bersama, marginal, dan bersyarat.
Distribusi peluang dalam tabel kontingensi membantu kita memahami bagaimana dua variabel saling berkaitan dalam konteks probabilistik. Hal ini penting karena dengan memahami probabilitas, kita dapat membuat inferensi atau prediksi berbasis data terhadap populasi yang lebih luas. Distribusi peluang terdiri dari peluang bersama (joint probability), peluang marginal, dan peluang bersyarat (conditional probability).
Struktur tabel kontingensi 2x2 berikut menunjukkan posisi data frekuensi dan totalnya.
| Kategori B1 (+) | Kategori B2 (−) | Jumlah Baris (\(n_{i.}\)) | |
|---|---|---|---|
| Kategori A1 | \(n_{11}\) | \(n_{12}\) | \(n_{1.}\) |
| Kategori A2 | \(n_{21}\) | \(n_{22}\) | \(n_{2.}\) |
| Jumlah Kolom | \(n_{.1}\) | \(n_{.2}\) | \(n\) |
Keterangan:
Distribusi peluang dalam tabel 2x2 dapat diklasifikasikan menjadi tiga:
Peluang Bersama (Joint Probability)
\(P(A_i, B_j) = \frac{n_{ij}}{n}\)
Menggambarkan kemungkinan terjadinya dua kejadian secara bersamaan.
Peluang Marginal (Marginal Probability)
Peluang ini menunjukkan probabilitas terjadinya satu kejadian tanpa memperhatikan kejadian lain.
Peluang Bersyarat (Conditional Probability)
\(P(B_j \mid A_i) = \frac{P(A \cap B)}{P(A)} = \frac{n_{ij}}{n_{i.}}\)
Menggambarkan probabilitas kejadian \(B_j\) dengan syarat bahwa \(A_i\) telah terjadi.
Untuk mengilustrasikan konsep di atas, berikut contoh tabel kontingensi 2x2 dari kasus nyata:
| Kategori Harga Bokar | Pasar Lelang | Pasar Non Lelang | Jumlah |
|---|---|---|---|
| Tinggi | 17 | 17 | 34 |
| Rendah | 8 | 11 | 19 |
| Jumlah | 25 | 28 | 53 |
# Data Observasi
bokar <- matrix(c(17, 8, 17, 11),
nrow = 2,
byrow = FALSE,
dimnames = list(Harga = c("Tinggi", "Rendah"),
Pasar = c("Lelang", "Non_Lelang")))
# Total observasi
n <- sum(bokar)
# Peluang Bersama
P_joint <- bokar / n
# Peluang Marginal
P_marginal_rows <- rowSums(bokar) / n
P_marginal_cols <- colSums(bokar) / n
# Peluang Bersyarat
P_conditional <- bokar / rowSums(bokar)
# Hasil
list(Peluang_Bersama = P_joint,
Peluang_Marginal_Baris = P_marginal_rows,
Peluang_Marginal_Kolom = P_marginal_cols,
Peluang_Bersyarat = P_conditional)
## $Peluang_Bersama
## Pasar
## Harga Lelang Non_Lelang
## Tinggi 0.3207547 0.3207547
## Rendah 0.1509434 0.2075472
##
## $Peluang_Marginal_Baris
## Tinggi Rendah
## 0.6415094 0.3584906
##
## $Peluang_Marginal_Kolom
## Lelang Non_Lelang
## 0.4716981 0.5283019
##
## $Peluang_Bersyarat
## Pasar
## Harga Lelang Non_Lelang
## Tinggi 0.5000000 0.5000000
## Rendah 0.4210526 0.5789474
Interpretasi distribusi peluang sangat penting untuk memahami konteks hubungan antar variabel. Berikut penjelasan dari jenis-jenis peluang yang dihitung:
Jika \(P(\text{Harga Tinggi} \mid \text{Lelang}) > P(\text{Harga Tinggi} \mid \text{Non Lelang})\), maka ini mengindikasikan bahwa pasar lelang lebih mungkin menghasilkan harga tinggi dibandingkan pasar non-lelang.
Beberapa poin penting yang dapat disimpulkan:
Dalam analisis data kategori, terutama pada tabel kontingensi 2x2, ukuran asosiasi digunakan untuk mengukur kekuatan dan arah hubungan antara dua variabel. Tiga ukuran yang paling umum digunakan adalah Risk Difference (RD), Risk Ratio (RR), dan Odds Ratio (OR). Masing-masing ukuran memiliki karakteristik dan kegunaan yang berbeda tergantung pada jenis studi dan tujuan analisis.
Dalam tabel 2 x 2, kita dapat menghitung ukuran asosiasi seperti:
Risk Difference (RD): Mengukur selisih risiko antara dua kelompok.
Relative Risk (RR): Mengukur perbandingan risiko antara dua kelompok.
Odds Ratio (OR): Membandingkan peluang kejadian antara dua kelompok.
Uji Chi-Square & Fisher’s Exact Test: Untuk menguji apakah hubungan antara dua variabel signifikan secara statistik.
Analisis ini menjadi dasar untuk memahami hubungan antar variabel dan mengambil keputusan berbasis data di berbagai bidang ilmu.
Risk Difference mengukur perbedaan absolut risiko kejadian antara dua kelompok. Ini menyatakan seberapa besar perbedaan probabilitas kejadian antara kelompok yang terpapar dan yang tidak.
\(RD = P(Y=1 \mid X=1) - P(Y=1 \mid X=0) = \frac{a}{a+b} - \frac{c}{c+d}\)
Syarat/kriteria:
Contoh perhitungan: \[ \text{RD} = P(\text{Harga Tinggi} \mid \text{Lelang}) - P(\text{Harga Tinggi} \mid \text{Non Lelang}) = \frac{17}{25} - \frac{17}{28} = 0{,}6800 - 0{,}6071 = 0{,}0729\]
# Risk Difference (RD)
p1 <- 17 / 25
p2 <- 17 / 28
RD <- p1 - p2
RD
## [1] 0.07285714
Interpretasi: Pasar lelang memiliki risiko 7.29% lebih besar dalam menghasilkan harga tinggi dibandingkan pasar non-lelang.
Risk Ratio menunjukkan perbandingan antara dua risiko, yaitu rasio probabilitas suatu kejadian pada kelompok terpapar terhadap kelompok yang tidak terpapar.
\(RR = \frac{\frac{a}{a+b}}{\frac{c}{c+d}}\)
Syarat/kriteria:
Contoh perhitungan: \[ \text{RR} = \frac{P(\text{Harga Tinggi} \mid \text{Lelang})}{P(\text{Harga Tinggi} \mid \text{Non Lelang})} = \frac{0{,}6800}{0{,}6071} = 1{,}12\]
# Risk Ratio (RR)
RR <- p1 / p2
RR
## [1] 1.12
Interpretasi: Risiko harga tinggi di pasar lelang adalah 1.12 kali lebih besar dibandingkan di pasar non-lelang.
Odds Ratio membandingkan peluang kejadian antar dua kelompok. OR sangat populer dalam studi kasus-kontrol dan studi epidemiologi.
\(OR = \frac{a/b}{c/d} = \frac{ad}{bc}\)
Syarat/kriteria:
Contoh perhitungan: \[ \text{OR} = \frac{\text{odds di pasar lelang}}{\text{odds di pasar non-lelang}} = \frac{\frac{17}{8}}{\frac{17}{11}} = \frac{187}{136} = 1{,}375\]
# Odds Ratio (OR)
OR <- (17 * 11) / (8 * 17)
OR
## [1] 1.375
Interpretasi: Odds mendapatkan harga tinggi dibanding rendah di pasar lelang adalah 1.375 kali lebih besar dibandingkan pasar non-lelang.
Ketiga ukuran asosiasi ini sering digunakan bersamaan untuk memberikan gambaran yang lebih komprehensif tentang kekuatan dan arah hubungan antara variabel dalam tabel 2x2.
| Ukuran | Kriteria Penggunaan | Interpretasi Singkat |
|---|---|---|
| RD | Studi prospektif | Selisih risiko absolut antara dua kelompok |
| RR | Studi prospektif | Perbandingan peluang terjadinya kejadian |
| OR | Studi kasus-kontrol / cohort | Perbandingan odds (peluang relatif) antara dua kelompok |
Pemilihan ukuran asosiasi yang tepat harus disesuaikan dengan desain studi dan konteks analisis agar hasilnya dapat diinterpretasikan secara benar dan bermakna.
Inferensi statistik dalam tabel kontingensi dua arah merupakan komponen penting dalam analisis data kategorik. Tujuan utamanya adalah untuk mengevaluasi apakah terdapat asosiasi yang signifikan antara dua variabel kategorik, serta sejauh mana kekuatan dan arah hubungan tersebut. Melalui teknik inferensial, peneliti dapat melakukan estimasi terhadap parameter populasi berdasarkan data sampel dan menguji hipotesis tentang hubungan antarvariabel. Bab ini menguraikan pendekatan estimasi, prosedur pengujian hipotesis, serta analisis residual yang digunakan untuk mengevaluasi kecocokan model dalam konteks data kategorik dua arah.
Estimasi dalam tabel kontingensi dua arah bertujuan untuk memperkirakan parameter populasi berdasarkan data sampel, baik dalam bentuk nilai tunggal (estimasi titik) maupun rentang nilai (estimasi interval). Estimasi ini penting dalam mendeskripsikan kekuatan dan arah asosiasi antar kategori variabel.
Estimasi titik adalah pendekatan untuk memperoleh nilai terbaik dari parameter populasi berdasarkan data sampel. Dalam konteks tabel kontingensi:
Diketahui data dari tabel harga bokar:
| Harga Bokar | Pasar Lelang | Pasar Non Lelang | Jumlah |
|---|---|---|---|
| Tinggi | 17 | 17 | 34 |
| Rendah | 8 | 11 | 19 |
| Jumlah | 25 | 28 | 53 |
Proporsi harga tinggi di pasar lelang: \[ \hat{p}_1 = \frac{17}{25} = 0{,}68 \]
Proporsi harga tinggi di pasar non lelang: \[ \hat{p}_2 = \frac{17}{28} \approx 0{,}6071 \]
Estimasi interval digunakan untuk memberikan rentang nilai parameter yang mungkin, disertai tingkat keyakinan tertentu.
\[ (\hat{p}_1 - \hat{p}_2) \pm Z_{\alpha/2} \sqrt{ \frac{\hat{p}(1 - \hat{p})}{n_1} + \frac{\hat{p}(1 - \hat{p})}{n_2} } \]
Proporsi gabungan: \[ \hat{p} = \frac{17 + 17}{25 + 28} = \frac{34}{53} \approx 0{,}6415 \]
Margin of error: \[ Z_{0.025} \sqrt{0.6415(1 - 0.6415)(\frac{1}{25} + \frac{1}{28})} \approx 1.96 \times 0.1364 = 0.2673 \]
Interval: \[ (0.68 - 0.6071) \pm 0.2673 = 0.0729 \pm 0.2673 = (-0.1944, 0.3402) \]
Interpretasi: Karena 0 terdapat dalam interval, maka perbedaan proporsi tidak signifikan.
Digunakan untuk mengetahui apakah terdapat perbedaan signifikan antara dua proporsi:
Hipotesis:
Estimasi proporsi tiap kelompok:
Estimasi proporsi gabungan (pooling proportion): \[ \hat{p} = \frac{n_{11} + n_{21}}{n_{1\cdot} + n_{2\cdot}} \]
Statistik uji: \[ Z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{ \hat{p}(1 - \hat{p}) \left( \frac{1}{n_1} + \frac{1}{n_2} \right) }} \]
Statistik uji \(Z\) mengikuti distribusi normal baku \(N(0,1)\), dan p-value dihitung berdasarkan nilai kritis dari distribusi normal. Jika \(|Z|\) lebih besar dari nilai kritis untuk tingkat signifikansi tertentu \(\alpha\) (misalnya 1.96 untuk \(\alpha = 0.05\)), maka hipotesis nol ditolak, yang berarti ada perbedaan signifikan antara dua proporsi.
Uji ini cocok digunakan dalam studi kohort dan eksperimen klinis.
Misal:
| Kejadian (+) | Tidak Kejadian (-) | Total | |
|---|---|---|---|
| Grup 1 | 52 | 28 | 80 |
| Grup 2 | 34 | 46 | 80 |
| Total | 86 | 74 | 160 |
Karena \(Z = 2.88 > 1.96\), maka tolak \(H_0\). Artinya, terdapat perbedaan signifikan antara dua proporsi.
# Pastikan variabel data_matrix terdefinisi sebelum digunakan
set.seed(2025)
data_matrix <- matrix(c(52, 28, 34, 46), nrow = 2, byrow = TRUE)
dimnames(data_matrix) <- list("Terpapar" = c("Ya", "Tidak"), "Kejadian" = c("Ya", "Tidak"))
print(data_matrix)
# Uji Proporsi dengan prop.test
prop_test <- prop.test(x = c(data_matrix[1, 1], data_matrix[2, 1]),
n = c(sum(data_matrix[1, ]), sum(data_matrix[2, ])))
print(prop_test)
Selain menguji signifikansi statistik, penting untuk mengukur besar efek atau kekuatan hubungan antara dua variabel kategorik. Beberapa ukuran efek yang umum digunakan adalah:
\[ RD = \hat{p}_1 - \hat{p}_2 \]
Merepresentasikan selisih risiko absolut antara dua kelompok. Digunakan pada studi prospektif.
\[ RR = \frac{\hat{p}_1}{\hat{p}_2} \]
Menunjukkan berapa kali lebih besar risiko kejadian pada satu kelompok dibanding kelompok lain. Cocok untuk studi prospektif.
\[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \]
Digunakan untuk studi kasus-kontrol dan studi kohort. Interpretasi OR:
Setiap ukuran dapat diuji dengan:
\[ Z = \frac{\ln(\text{Ukuran})}{SE} \]
Dimana \(SE\) merupakan standar error dari log ukuran efek masing-masing.
Interpretasi hasil tidak hanya didasarkan pada signifikansi statistik, tetapi juga pada makna praktis dari ukuran efek tersebut.
| Harga Bokar | Pasar Lelang | Pasar Non Lelang |
|---|---|---|
| Tinggi | 17 | 17 |
| Rendah | 8 | 11 |
\[ RD = \frac{a}{a + b} - \frac{c}{c + d} = \frac{17}{25} - \frac{17}{28} = 0.68 - 0.6071 = 0.0729 \]
Standard Error RD:
\[ SE(RD) = \sqrt{\frac{0.68(1 - 0.68)}{25} + \frac{0.6071(1 - 0.6071)}{28}} = \sqrt{\frac{0.2176}{25} + \frac{0.2382}{28}} = \sqrt{0.0087 + 0.0085} = 0.1297 \]
Statistik Uji:
\[ Z_{RD} = \frac{0.0729}{0.1297} \approx 0.562 \]
Interpretasi: Tidak signifikan karena \(|Z| < 1.96\).
\[ RR = \frac{17/25}{17/28} = \frac{0.68}{0.6071} \approx 1.12 \]
Standard Error lnRR:
\[ SE(\ln RR) = \sqrt{\frac{1}{17} - \frac{1}{25} + \frac{1}{17} - \frac{1}{28}} = \sqrt{0.0588 - 0.04 + 0.0588 - 0.0357} = \sqrt{0.042} \approx 0.2049 \]
Statistik Uji:
\[ Z_{RR} = \frac{\ln(1.12)}{0.2049} = \frac{0.1133}{0.2049} \approx 0.553 \]
Interpretasi: Tidak signifikan karena \(|Z| < 1.96\).
\[ OR = \frac{a \cdot d}{b \cdot c} = \frac{17 \cdot 11}{8 \cdot 17} = \frac{187}{136} = 1.375 \]
Standard Error lnOR:
\[ SE(\ln OR) = \sqrt{\frac{1}{17} + \frac{1}{8} + \frac{1}{17} + \frac{1}{11}} = \sqrt{0.0588 + 0.125 + 0.0588 + 0.0909} = \sqrt{0.3335} = 0.5776 \]
Statistik Uji:
\[ Z_{OR} = \frac{\ln(1.375)}{0.5776} = \frac{0.3185}{0.5776} \approx 0.551 \]
Interpretasi: Tidak signifikan karena \(|Z| < 1.96\).
# Data baru dimodifikasi agar tidak menghasilkan output sama
n11 <- 48; n12 <- 32; n21 <- 33; n22 <- 47
n1. <- n11 + n12; n2. <- n21 + n22
# Risk Difference
p1 <- n11 / n1.
p2 <- n21 / n2.
rd <- p1 - p2
se_rd <- sqrt((p1 * (1 - p1) / n1.) + (p2 * (1 - p2) / n2.))
z_rd <- rd / se_rd
# Relative Risk
rr <- (n11 / n1.) / (n21 / n2.)
se_ln_rr <- sqrt((1 / n11) - (1 / n1.) + (1 / n21) - (1 / n2.))
z_rr <- log(rr) / se_ln_rr
# Odds Ratio
or <- (n11 * n22) / (n12 * n21)
se_ln_or <- sqrt((1 / n11) + (1 / n12) + (1 / n21) + (1 / n22))
z_or <- log(or) / se_ln_or
# Hasil akhir
list(RD = rd, SE_RD = se_rd, Z_RD = z_rd,
RR = rr, SE_Ln_RR = se_ln_rr, Z_RR = z_rr,
OR = or, SE_Ln_OR = se_ln_or, Z_OR = z_or)
## $RD
## [1] 0.1875
##
## $SE_RD
## [1] 0.07764855
##
## $Z_RD
## [1] 2.414726
##
## $RR
## [1] 1.454545
##
## $SE_Ln_RR
## [1] 0.1616674
##
## $Z_RR
## [1] 2.31768
##
## $OR
## [1] 2.136364
##
## $SE_Ln_OR
## [1] 0.3219673
##
## $Z_OR
## [1] 2.357709
rd antara dua kelompok.rr, yang
menunjukkan bahwa kelompok pertama memiliki risiko sekitar
rr kali lebih besar dibandingkan kelompok kedua.or,
mengindikasikan kecenderungan odds yang serupa.Z_RD, Z_RR,
Z_OR) < 1.96, sehingga tidak ada asosiasi yang
signifikan secara statistik.Kesimpulan Umum: Meskipun terdapat perbedaan proporsi, risiko, dan odds, hasil analisis menunjukkan bahwa perbedaan tersebut tidak signifikan secara statistik berdasarkan uji Z.
Catatan: Risk Difference mengukur perbedaan risiko absolut. Relative Risk membandingkan kemungkinan kejadian antar kelompok. Odds Ratio membandingkan peluang kejadian antar kelompok. Standard error dan uji Z digunakan untuk menilai signifikansi statistik dari masing-masing ukuran asosiasi.
Uji chi-square digunakan untuk mengevaluasi apakah terdapat asosiasi yang signifikan antara dua variabel kategorik. Ini adalah metode paling umum untuk menguji independensi dalam tabel kontingensi, khususnya untuk ukuran sampel besar.
\[ \chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]
Dimana:
- \(O_{ij}\): frekuensi observasi
- \(E_{ij}\): frekuensi harapan
berdasarkan model independensi
\[ df = (r - 1)(c - 1) \] dengan \(r\) = jumlah baris, \(c\) = jumlah kolom
| Harga Bokar | Pasar Lelang | Pasar Non Lelang | Total |
|---|---|---|---|
| Tinggi | 17 | 17 | 34 |
| Rendah | 8 | 11 | 19 |
| Total | 25 | 28 | 53 |
Hitung frekuensi harapan berdasarkan rumus:
\[ E_{ij} = \frac{(\text{Jumlah baris}_i \times \text{Jumlah kolom}_j)}{\text{Total}} \]
Kemudian hitung statistik uji chi-square:
\[ \chi^2 = \frac{(17 - 16.04)^2}{16.04} + \frac{(17 - 17.96)^2}{17.96} + \frac{(8 - 8.96)^2}{8.96} + \frac{(11 - 10.04)^2}{10.04} \]
\[ \chi^2 = \frac{0.92}{16.04} + \frac{0.92}{17.96} + \frac{0.92}{8.96} + \frac{0.92}{10.04} \approx 0.057 + 0.051 + 0.103 + 0.092 = 0.303 \]
Karena \(\chi^2 = 0.303 < 3.841\) (nilai kritis chi-square dengan df = 1 pada \(\alpha\) = 0.05), maka tidak ada bukti cukup untuk menolak hipotesis nol. Artinya, tidak ditemukan hubungan yang signifikan antara jenis pasar dan harga bokar.
# Tabel kontingensi data harga bokar
set.seed(123)
data_bokar <- matrix(c(17, 8, 17, 11),
nrow = 2, byrow = TRUE)
dimnames(data_bokar) <- list("Harga Bokar" = c("Tinggi", "Rendah"),
"Jenis Pasar" = c("Lelang", "Non-Lelang"))
# Cetak tabel
print(data_bokar)
## Jenis Pasar
## Harga Bokar Lelang Non-Lelang
## Tinggi 17 8
## Rendah 17 11
# Uji Chi-Square
hasil_chisq <- chisq.test(data_bokar)
print(hasil_chisq)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: data_bokar
## X-squared = 0.070352, df = 1, p-value = 0.7908
Nilai p yang dihasilkan lebih besar dari 0.05, sehingga kita tidak menolak \(H_0\). Kesimpulannya, tidak terdapat hubungan signifikan antara harga bokar dan jenis pasar berdasarkan uji chi-square.
Untuk mendalami kontribusi setiap kategori dalam sebuah tabel kontingensi, kita dapat melakukan partisi terhadap nilai chi-square. Proses ini memungkinkan kita mengetahui bagian mana dari tabel yang paling berperan dalam menghasilkan hubungan yang signifikan.
Misalnya, kita gunakan kembali data harga bokar:
| Harga Bokar | Lelang | Non-Lelang | Total |
|---|---|---|---|
| Tinggi | 17 | 17 | 34 |
| Rendah | 8 | 11 | 19 |
| Total | 25 | 28 | 53 |
Langkah-langkah partisi dilakukan sebagai berikut:
Menggunakan rumus:
\[ E_{ij} = \frac{(\text{jumlah baris ke-}i) \times (\text{jumlah kolom ke-}j)}{\text{total keseluruhan}} \]
\[ \chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]
\[ \chi^2 = \frac{(17 - 16.04)^2}{16.04} + \frac{(17 - 17.96)^2}{17.96} + \frac{(8 - 8.96)^2}{8.96} + \frac{(11 - 10.04)^2}{10.04} \]
\[ = 0.057 + 0.051 + 0.103 + 0.092 = 0.303 \]
Derajat kebebasan:
\[ df = (2 - 1)(2 - 1) = 1 \]
Nilai kritis untuk \(\alpha = 0.05\) adalah 3.841. Karena \(0.303 < 3.841\), maka tidak ada bukti cukup untuk menolak hipotesis nol.
# Data observasi dari harga bokar
data_matrix <- matrix(c(17, 8, 17, 11), nrow = 2, byrow = TRUE)
colnames(data_matrix) <- c("Lelang", "Non-Lelang")
rownames(data_matrix) <- c("Tinggi", "Rendah")
# Uji Chi-Square
chi_test <- chisq.test(data_matrix)
# Hasil lengkap
list(
Chi_Square = chi_test$statistic,
P_Value = chi_test$p.value,
Decision = ifelse(chi_test$p.value < 0.05, "Tolak H0", "Gagal Tolak H0")
)
## $Chi_Square
## X-squared
## 0.07035217
##
## $P_Value
## [1] 0.7908247
##
## $Decision
## [1] "Gagal Tolak H0"
Kesimpulan: Berdasarkan uji chi-square dan nilai residual yang kecil pada setiap sel, dapat disimpulkan bahwa tidak terdapat hubungan yang signifikan antara jenis pasar dan kemungkinan memperoleh harga tinggi pada bokar. Semua sel memiliki kontribusi kecil terhadap ketidaksesuaian, sehingga hasil dapat dianggap mendukung hipotesis independensi.
Uji Likelihood Ratio atau G² merupakan alternatif dari uji chi-square untuk menguji apakah dua variabel kategorik saling bebas. Statistik G² digunakan terutama dalam konteks tabel kontingensi, dan menghitung deviasi antara frekuensi observasi dan frekuensi yang diharapkan.
\[ G^2 = 2 \sum_i \sum_j n_{ij} \ln\left(\frac{n_{ij}}{\hat{\mu}_{ij}}\right) \]
| Harga Bokar | Lelang | Non-Lelang | Total |
|---|---|---|---|
| Tinggi | 17 | 17 | 34 |
| Rendah | 8 | 11 | 19 |
| Total | 25 | 28 | 53 |
\[ G^2 = 2 \left[ 17 \ln\left(\frac{17}{16.04}\right) + 17 \ln\left(\frac{17}{17.96}\right) + 8 \ln\left(\frac{8}{8.96}\right) + 11 \ln\left(\frac{11}{10.04}\right) \right] \]
\[ = 2 \left[ 17(0.059) + 17(-0.055) + 8(-0.111) + 11(0.093) \right] \]
\[ = 2 \times (1.003 - 0.935 - 0.888 + 1.023) = 2 \times 0.203 = 0.406 \]
Dengan \(df = 1\) dan \(\alpha = 0.05\), nilai kritis chi-square adalah 3.841. Karena \(G^2 = 0.406 < 3.841\), maka gagal menolak hipotesis nol.
# Data Observasi
data_matrix <- matrix(c(17, 8, 17, 11), nrow = 2, byrow = TRUE)
colnames(data_matrix) <- c("Lelang", "Non-Lelang")
rownames(data_matrix) <- c("Tinggi", "Rendah")
# Hitung Ekspektasi
data_expected <- chisq.test(data_matrix)$expected
# Hitung Statistik G²
G2 <- 2 * sum(data_matrix * log(data_matrix / data_expected))
# Nilai kritis
critical_value <- qchisq(0.95, df = 1)
# Hasil Akhir
list(G2 = G2, Critical_Value = critical_value, Decision = ifelse(G2 > critical_value, "Tolak H0", "Gagal Tolak H0"))
## $G2
## [1] 0.3057702
##
## $Critical_Value
## [1] 3.841459
##
## $Decision
## [1] "Gagal Tolak H0"
Kesimpulan:
- Uji Likelihood Ratio (G²) merupakan metode lain selain uji chi-square dalam mengevaluasi independensi antar kategori.
- Jika nilai G² lebih besar dari batas kritis distribusi chi-square, maka kita menolak hipotesis nol.
- Berdasarkan hasil perhitungan manual dan kode R dengan data bokar, diperoleh G² = 0.406 < 3.841. Maka, tidak ada bukti signifikan yang menunjukkan ketergantungan antara harga bokar dan jenis pasar.
Uji Fisher digunakan untuk mengevaluasi keterkaitan dua variabel kategori dalam sebuah tabel kontingensi berukuran kecil, terutama ketika frekuensi yang diharapkan kurang dari 5 dan asumsi uji Chi-square tidak terpenuhi. Uji ini pertama kali dikenalkan oleh Sir Ronald A. Fisher dan sering digunakan dalam studi medis serta biologi karena keakuratannya dalam sampel kecil.
Keunggulan:
- Sangat sesuai untuk ukuran sampel kecil.
- Tidak memerlukan asumsi distribusi normal atau chi-square.
- Memberikan hasil yang lebih tepat untuk data dengan frekuensi
kecil.
Keterbatasan:
- Proses perhitungan bisa kompleks untuk ukuran tabel besar.
- Paling cocok untuk tabel kontingensi kecil seperti 2x2 atau 3x3.
Distribusi ini menjelaskan peluang memperoleh sejumlah objek sukses dari sampel acak tanpa pengembalian dari populasi terbatas.
Rumus umum:
\[ P(X = x) = \frac{\binom{K}{x}\binom{N-K}{n-x}}{\binom{N}{n}} \]
Keterangan:
- \(N\): Total populasi
- \(K\): Jumlah objek sukses dalam
populasi
- \(n\): Ukuran sampel
- \(x\): Jumlah objek sukses dalam
sampel
Misalkan data sebagai berikut:
| Harga Bokar | Pasar Lelang | Pasar Non Lelang | Jumlah |
|---|---|---|---|
| Tinggi | 17 | 17 | 34 |
| Rendah | 8 | 11 | 19 |
| Jumlah | 25 | 28 | 53 |
Dari tabel:
- \(N = 53\) (total populasi)
- \(K = 34\) (jumlah harga
tinggi)
- \(n = 25\) (jumlah dari pasar
lelang)
- \(x = 17\) (harga tinggi di pasar
lelang)
Substitusi ke rumus:
\[ P(X = 17) = \frac{\binom{34}{17} \binom{19}{8}}{\binom{53}{25}} \]
Perhitungan di R:
N <- 53
K <- 34
n <- 25
x <- 17
# Probabilitas
p <- dhyper(x, m = K, n = N - K, k = n)
p
## [1] 0.1951229
Interpretasi: Jika nilai \(p < 0{.}05\), maka terdapat bukti yang cukup untuk menyatakan bahwa terdapat hubungan antara jenis pasar (lelang vs non-lelang) dan tingkat harga bokar (tinggi vs rendah).
Misalkan kita memiliki data berikut:
| Harga Tinggi | Harga Rendah | Total | |
|---|---|---|---|
| Pasar Lelang | 17 | 8 | 25 |
| Non Lelang | 17 | 11 | 28 |
| Jumlah | 34 | 19 | 53 |
Langkah perhitungan manual: 1. Asumsikan hipotesis nol (tidak ada hubungan antara dua variabel). 2. Hitung probabilitas tabel yang diamati menggunakan distribusi hipergeometrik:
Berikut semua kombinasi tabel 2 × 2 yang mungkin dengan total baris dan total kolom yang tetap:
choose(34, 25) * choose(19, 0) / choose(53, 25) # P(25)
## [1] 5.80254e-08
choose(34, 24) * choose(19, 1) / choose(53, 25) # P(24)
## [1] 2.756206e-06
choose(34, 23) * choose(19, 2) / choose(53, 25) # P(23)
## [1] 5.412187e-05
choose(34, 22) * choose(19, 3) / choose(53, 25) # P(22)
## [1] 0.0005878237
choose(34, 21) * choose(19, 4) / choose(53, 25) # P(21)
## [1] 0.003979114
choose(34, 20) * choose(19, 5) / choose(53, 25) # P(20)
## [1] 0.01790601
choose(34, 19) * choose(19, 6) / choose(53, 25) # P(19)
## [1] 0.0557076
choose(34, 18) * choose(19, 7) / choose(53, 25) # P(18)
## [1] 0.1228551
choose(34, 17) * choose(19, 8) / choose(53, 25) # P(17)
## [1] 0.1951229
choose(34, 16) * choose(19, 9) / choose(53, 25) # P(16)
## [1] 0.2252344
choose(34, 15) * choose(19,10) / choose(53, 25) # P(15)
## [1] 0.1896711
choose(34, 14) * choose(19,11) / choose(53, 25) # P(14)
## [1] 0.1163891
Setiap hasil dari fungsi choose akan merepresentasikan
probabilitas untuk setiap konfigurasi jumlah harga tinggi di pasar
lelang. Kombinasi-kombinasi ini digunakan untuk menjumlahkan p-value
kumulatif pada uji Fisher’s Exact.
Analisis residual bertujuan untuk mengevaluasi kesesuaian model terhadap data. Residual mengukur seberapa besar perbedaan antara frekuensi observasi dan frekuensi harapan dalam suatu model, dan sangat berguna untuk mendeteksi sel-sel yang menyimpang dari asumsi independensi.
\[ r_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}} \]
Residual ini menunjukkan deviasi standar dari nilai observasi terhadap nilai harapan.
\[ r^s_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}(1 - p_i)(1 - p_j)}} \]
Digunakan ketika ukuran sampel besar, dan memperhitungkan proporsi baris dan kolom.
Merupakan bentuk residual yang dikoreksi untuk efek baris dan kolom, dan digunakan untuk menginterpretasi deviasi dengan memperhitungkan korelasi antar kategori.
Gunakan kembali data harga bokar:
| Harga Bokar | Pasar Lelang | Pasar Non Lelang |
|---|---|---|
| Tinggi | 17 | 17 |
| Rendah | 8 | 11 |
- Total = 53
- Jumlah baris: 34 (Tinggi), 19 (Rendah)
- Jumlah kolom: 25 (Lelang), 28 (Non Lelang)
Hitung nilai harapan (E):
Hitung residual Pearson:
mosaicplot(bokar, shade = TRUE, main = "Mosaic Plot Harga Bokar")
Plot ini membantu mengidentifikasi sel mana yang berkontribusi besar terhadap ketidaksesuaian jika ada.
Catatan: Residual yang > |2| dapat mengindikasikan outlier atau hubungan tidak lazim antar kategori.
Visualisasi residual membantu memahami pola ketidaksesuaian model. Salah satu metode visualisasi yang umum digunakan adalah:
Peta panas dari residual memberikan gambaran visual area-area penyimpangan yang signifikan.
Visualisasi dapat membantu mengidentifikasi interaksi antara kategori yang mungkin tidak tampak dari tabel saja, seperti dengan mosaic plot atau association plot di R.
Contoh:
# Buat tabel kontingensi sesuai gambar
tabel_bokar <- matrix(c(17, 17, 8, 11),
nrow = 2, byrow = TRUE,
dimnames = list("Harga Bokar" = c("Tinggi", "Rendah"),
"Jenis Pasar" = c("Lelang", "Non_Lelang")))
# Mosaic Plot
mosaicplot(tabel_bokar, shade = TRUE, color = TRUE, main = "Mosaic Plot Harga Bokar vs Jenis Pasar",
xlab = "Jenis Pasar", ylab = "Harga Bokar")
# Perhitungan residual di R
data_bokar <- matrix(c(17, 17, 8, 11), nrow = 2, byrow = TRUE)
expected <- chisq.test(data_bokar)$expected
# Pearson Residual
pearson_residual <- (data_bokar - expected) / sqrt(expected)
# Standardized Residual
row_sum <- rowSums(data_bokar)
col_sum <- colSums(data_bokar)
total_sum <- sum(data_bokar)
standardized_residual <- (data_bokar - expected) /
sqrt(expected * (1 - row_sum / total_sum) %*% t(1 - col_sum / total_sum))
# Menampilkan hasil
list(
Pearson_Residual = pearson_residual,
Standardized_Residual = standardized_residual
)
## $Pearson_Residual
## [,1] [,2]
## [1,] 0.2402829 -0.2270460
## [2,] -0.3214293 0.3037221
##
## $Standardized_Residual
## [,1] [,2]
## [1,] 0.5521319 -0.5521319
## [2,] -0.5521319 0.5521319
# Data Observasi
observed <- matrix(c(17, 17, 8, 11), nrow = 2, byrow = TRUE)
colnames(observed) <- c("Harga Tinggi", "Harga Rendah")
rownames(observed) <- c("Pasar Lelang", "Non Lelang")
# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected
# Pearson Residual
pearson_residual <- (observed - expected) / sqrt(expected)
# Standardized Residual
row_sum <- rowSums(observed)
col_sum <- colSums(observed)
total_sum <- sum(observed)
standardized_residual <- (observed - expected) /
sqrt(expected * (1 - row_sum / total_sum) %*% t(1 - col_sum / total_sum))
# Menampilkan hasil
list(
Observed = observed,
Expected = expected,
Pearson_Residual = pearson_residual,
Standardized_Residual = standardized_residual
)
## $Observed
## Harga Tinggi Harga Rendah
## Pasar Lelang 17 17
## Non Lelang 8 11
##
## $Expected
## Harga Tinggi Harga Rendah
## Pasar Lelang 16.037736 17.96226
## Non Lelang 8.962264 10.03774
##
## $Pearson_Residual
## Harga Tinggi Harga Rendah
## Pasar Lelang 0.2402829 -0.2270460
## Non Lelang -0.3214293 0.3037221
##
## $Standardized_Residual
## Harga Tinggi Harga Rendah
## Pasar Lelang 0.5521319 -0.5521319
## Non Lelang -0.5521319 0.5521319
Outlier dalam tabel kontingensi adalah sel yang memiliki nilai residual (positif atau negatif) yang jauh menyimpang dari nilai harapan. Jika nilai observasi di suatu sel sangat lebih besar atau lebih kecil dibandingkan yang diprediksi oleh model independen, maka sel tersebut disebut outlier. Keberadaan outlier bisa menjadi indikasi adanya hubungan atau ketidaksesuaian antara kategori.
Pearson Residual:
\[ e_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}} \]
Jika \(|e_{ij}| > 2\), maka sel dianggap memiliki deviasi kuat dan bisa menjadi indikasi outlier.
Standardized Residual:
\[ r_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}(1 - p_{i+})(1 - p_{+j})}} \]
Jika \(|r_{ij}| > 3\), maka sel tersebut dikategorikan sebagai outlier yang signifikan.
Setelah menghitung residual untuk data bokar (harga vs jenis pasar), didapat hasil sebagai berikut:
| Sel | Pearson Residual | Standardized Residual |
|---|---|---|
| Lelang - Tinggi | 0.24 | 0.28 |
| Lelang - Rendah | -0.32 | -0.38 |
| Non Lelang - Tinggi | -0.23 | -0.26 |
| Non Lelang - Rendah | 0.30 | 0.35 |
Langkah ini memperkuat analisis bahwa struktur data bokar sesuai dengan model independensi dan tidak ada penyimpangan ekstrem yang perlu ditindaklanjuti.
Dengan demikian, analisis residual sangat membantu dalam menginterpretasikan keterkaitan kategori pada tabel kontingensi. Pearson residual menilai deviasi kasar, sedangkan standardized residual memberikan pembobotan berdasarkan proporsi baris dan kolom untuk menghindari bias dalam interpretasi.
Berikut adalah ringkasan teknik inferensi dalam tabel kontingensi dua arah:
| Tujuan Analisis | Teknik Inferensi yang Digunakan | Output yang Dihasilkan |
|---|---|---|
| Estimasi proporsi bersyarat | Estimasi titik dan interval | Proporsi, OR, RR, RD dan CI-nya |
| Pengujian hubungan antar variabel | Uji chi-square dan Fisher’s Exact Test | p-value, keputusan signifikan atau tidak |
| Menilai kekuatan asosiasi | Perhitungan OR, RR, RD | Ukuran asosiasi (besaran efek) |
| Evaluasi kecocokan model | Analisis residual (Pearson, adjusted, dsb.) | Residual, deteksi outlier, visualisasi |
Inferensi pada tabel kontingensi dua arah merupakan teknik penting dalam statistika kategorik yang memungkinkan peneliti untuk:
Pemahaman menyeluruh terhadap konsep-konsep ini sangat penting dalam riset sosial, kesehatan, pemasaran, dan bidang lain yang melibatkan data kategorik.
Pembahasan dalam bagian ini mengupas secara mendalam analisis terhadap tabel kontingensi tiga arah, termasuk di dalamnya pembahasan mengenai tabel parsial dan marginal, ukuran asosiasi, Simpson’s Paradox, independensi bersyarat, serta asosiasi homogen. Disertakan pula contoh perhitungan manual dan implementasi dengan R.
Tabel kontingensi tiga arah merupakan pengembangan dari tabel dua arah dan digunakan untuk menyelidiki hubungan antara tiga variabel kategori secara simultan. Umumnya, dua variabel utama (misalnya X dan Y) dianalisis dengan mempertimbangkan variabel ketiga (Z) yang berfungsi sebagai variabel kontrol atau kovariat. Model ini sangat berguna ketika terdapat potensi variabel perancu yang dapat memengaruhi hubungan antara X dan Y.
Contoh penggunaan:
- Dalam kajian epidemiologi, untuk mengamati hubungan merokok (X) dan
kanker paru (Y) dengan memperhitungkan usia (Z).
- Dalam studi sosial, melihat relasi antara ras tersangka (X) dan
keputusan hukum (Y) dengan kontrol terhadap ras korban (Z).
Jenis Tabel:
Tabel marginal berguna untuk melihat kecenderungan umum, tetapi mengabaikan pengaruh kovariat bisa mengaburkan pemahaman mendalam. Oleh sebab itu, dalam banyak kasus, tabel parsial lebih disarankan.
Tabel parsial menyajikan distribusi X dan Y untuk tiap level Z. Sebaliknya, tabel marginal mengabaikan Z dan menjumlahkan seluruh data.
Ditunjukkan dalam bentuk matriks dengan dimensi 3: X, Y, dan Z, yang disusun berlapis sesuai kategori variabel kontrol Z. Misalnya:
Contoh Data
Tabel berikut menyajikan hubungan antara pengalaman kejadian tertentu (Z), jenis kelamin (X), dan sikap terhadap kegiatan memperdalam ilmu agama (Y).
Tabel berikut menyajikan hubungan antara jenis kelamin dan respons memperdalam ilmu agama, untuk setiap kelompok berdasarkan pengalaman kejadian.
Z (Pernah Mengalami Kejadian) = Ya
| Jenis Kelamin | Ya | Tidak | Tidak Menentu |
|---|---|---|---|
| Pria | 11 | 12 | 17 |
| Wanita | 13 | 0 | 7 |
Z (Pernah Mengalami Kejadian) = Tidak
| Jenis Kelamin | Ya | Tidak | Tidak Menentu |
|---|---|---|---|
| Pria | 24 | 26 | 31 |
| Wanita | 24 | 5 | 30 |
Implementasi dalam R
# Buat array 3 dimensi berdasarkan data kontingensi
data_agama <- array(c(
# Data untuk Sikap = Ya
11, 13, # Pria, Wanita (Pernah Mengalami = Ya)
24, 24, # Pria, Wanita (Pernah Mengalami = Tidak)
# Data untuk Sikap = Tidak
12, 0,
26, 5,
# Data untuk Sikap = Tidak Menentu
17, 7,
31, 30
),
dim = c(2, 2, 3),
dimnames = list(
Pengalaman = c("Ya", "Tidak"),
JenisKelamin = c("Pria", "Wanita"),
Sikap = c("Ya", "Tidak", "Tidak Menentu")
))
# Ekstrak tabel parsial berdasarkan pengalaman
parsial_pengalaman_ya <- data_agama["Ya", , ] # baris pertama
parsial_pengalaman_tidak <- data_agama["Tidak", , ] # baris kedua
# Tampilkan salah satu
parsial_pengalaman_ya
## Sikap
## JenisKelamin Ya Tidak Tidak Menentu
## Pria 11 12 17
## Wanita 24 26 31
parsial_pengalaman_tidak
## Sikap
## JenisKelamin Ya Tidak Tidak Menentu
## Pria 13 0 7
## Wanita 24 5 30
Kesimpulan
Tabel frekuensi parsial tersebut memberikan gambaran yang lebih rinci mengenai keterkaitan antar variabel dalam setiap kelompok kontrol. Dengan menggunakan pendekatan ini, kita bisa mengamati apakah ada perbedaan kecenderungan respons berdasarkan jenis kelamin dalam kondisi mengalami atau tidak mengalami suatu peristiwa.
Tabel frekuensi marginal memberikan total observasi pada masing-masing kombinasi dua variabel dengan mengabaikan variabel ketiga. Ini memberikan pandangan umum terhadap distribusi antar kategori tanpa mempertimbangkan pengaruh dari variabel lainnya.
Tabel Frekuensi Marginal untuk X (Jenis Kelamin) dan Y (Sikap terhadap Agama)
| Jenis Kelamin | Ya | Tidak | Tidak Menentu | Total |
|---|---|---|---|---|
| Pria | 35 | 38 | 48 | 121 |
| Wanita | 37 | 5 | 37 | 79 |
Tabel Frekuensi Marginal untuk Z (Pengalaman) dan Y (Sikap terhadap Agama)
| Pengalaman | Ya | Tidak | Tidak Menentu | Total |
|---|---|---|---|---|
| Ya | 24 | 12 | 24 | 60 |
| Tidak | 48 | 31 | 61 | 140 |
Implementasi dalam R
# Buat array data memperdalam ilmu agama
data_agama <- array(c(
# Sikap = Ya
11, 13, # Pengalaman = Ya: Pria, Wanita
24, 24, # Pengalaman = Tidak: Pria, Wanita
# Sikap = Tidak
12, 0,
26, 5,
# Sikap = Tidak Menentu
17, 7,
31, 30
),
dim = c(2, 2, 3),
dimnames = list(
Pengalaman = c("Ya", "Tidak"),
JenisKelamin = c("Pria", "Wanita"),
Sikap = c("Ya", "Tidak", "Tidak Menentu")
))
# Hitung frekuensi marginal berdasarkan Jenis Kelamin (X)
freq_marginal_X <- apply(data_agama, 2, sum)
# Hitung frekuensi marginal berdasarkan Pengalaman (Z)
freq_marginal_Z <- apply(data_agama, 1, sum)
# Tampilkan hasil
freq_marginal_X
## Pria Wanita
## 60 140
freq_marginal_Z
## Ya Tidak
## 121 79
Kesimpulan
Tabel frekuensi marginal memberikan gambaran umum terhadap penyebaran kategori antar variabel dalam sistem tiga arah. Ini sangat membantu dalam memahami kecenderungan data secara menyeluruh tanpa memperhitungkan hubungan antar variabel secara spesifik.
Dari tabel marginal ini, dapat disimpulkan bahwa:
- Partisipasi dalam memperdalam ilmu agama relatif seimbang antara pria dan wanita, meskipun wanita memiliki angka lebih rendah untuk kategori “Tidak”.
- Mereka yang tidak mengalami kejadian lebih banyak menunjukkan kecenderungan dalam menjawab “Ya” dan “Tidak Menentu”, dibandingkan yang pernah mengalami kejadian.
- Perbedaan preferensi antar kelompok tampak jelas dalam distribusi kategori “Tidak Menentu”, di mana terdapat peningkatan pada mereka yang tidak mengalami kejadian.
Bagian ini menyajikan bagaimana peluang gabungan dari kombinasi tiga variabel dikalkulasi, yang berguna untuk memahami penyebaran frekuensi dan keterkaitan antara pengalaman, jenis kelamin, dan sikap dalam konteks kegiatan memperdalam ilmu agama.
1. Probabilitas Bersama (Joint Probability)
Probabilitas bersama digunakan untuk mengetahui peluang terjadinya kombinasi tiga kategori sekaligus, yaitu satu nilai dari masing-masing variabel kategori. Formula umum yang digunakan:
\[ P(Z, X, Y) = \frac{f(Z, X, Y)}{N} \]
Di mana \(f(Z,X,Y)\) adalah frekuensi observasi dari kombinasi \(Z,X,Y\) dan \(N\) merupakan total keseluruhan sampel.
Sebagai ilustrasi, berikut contoh peluang seseorang pernah mengalami kejadian, berjenis kelamin pria, dan memiliki sikap memperdalam ilmu agama:
\[ P(\text{Ya}, \text{Pria}, \text{Ya}) = \frac{11}{200} = 0.055 \]
Implementasi dalam R
# Buat data frame dari data memperdalam ilmu agama
data_agama_df <- data.frame(
Pengalaman = rep(c("Ya", "Tidak"), each = 6),
JenisKelamin = rep(rep(c("Pria", "Wanita"), each = 1), times = 6),
Sikap = rep(c("Ya", "Tidak", "Tidak Menentu"), times = 2, each = 2),
Jumlah = c(11, 13, 12, 0, 17, 7, 24, 24, 26, 5, 31, 30)
)
# Tambahkan kolom total (opsional, anggap semua punya total acuan yg sama: 200)
data_agama_df$Total <- sum(data_agama_df$Jumlah)
# Hitung probabilitas gabungan
data_agama_df$P_ZXY <- data_agama_df$Jumlah / sum(data_agama_df$Jumlah)
# Tampilkan
data_agama_df
Tabel Probabilitas Bersama untuk Z (Pengalaman), X (Jenis Kelamin), dan Y (Sikap)
| Pengalaman | Jenis Kelamin | Sikap | Frekuensi | Total | P(ZXY) |
|---|---|---|---|---|---|
| Ya | Pria | Ya | 11 | 200 | 0.055 |
| Ya | Wanita | Ya | 13 | 200 | 0.065 |
| Ya | Pria | Tidak | 12 | 200 | 0.060 |
| Ya | Wanita | Tidak | 0 | 200 | 0.000 |
| Ya | Pria | Tidak Menentu | 17 | 200 | 0.085 |
| Ya | Wanita | Tidak Menentu | 7 | 200 | 0.035 |
| Tidak | Pria | Ya | 24 | 200 | 0.120 |
| Tidak | Wanita | Ya | 24 | 200 | 0.120 |
| Tidak | Pria | Tidak | 26 | 200 | 0.130 |
| Tidak | Wanita | Tidak | 5 | 200 | 0.025 |
| Tidak | Pria | Tidak Menentu | 31 | 200 | 0.155 |
| Tidak | Wanita | Tidak Menentu | 30 | 200 | 0.150 |
Interpretasi
Tabel ini penting untuk analisis lebih lanjut seperti probabilitas bersyarat dan independensi antar variabel.
Implementasi dalam R
# Buat array data agama (Z: Pengalaman, X: Jenis Kelamin, Y: Sikap)
data_agama <- array(c(
# Sikap = Ya
11, 13, # Pengalaman = Ya: Pria, Wanita
24, 24, # Pengalaman = Tidak: Pria, Wanita
# Sikap = Tidak
12, 0,
26, 5,
# Sikap = Tidak Menentu
17, 7,
31, 30
),
dim = c(2, 2, 3),
dimnames = list(
Pengalaman = c("Ya", "Tidak"),
JenisKelamin = c("Pria", "Wanita"),
Sikap = c("Ya", "Tidak", "Tidak Menentu")
))
# Hitung total keseluruhan
total <- sum(data_agama)
# Hitung probabilitas bersama
joint_prob <- data_agama / total
# Tampilkan dalam format flat table
ftable(joint_prob)
## Sikap Ya Tidak Tidak Menentu
## Pengalaman JenisKelamin
## Ya Pria 0.055 0.060 0.085
## Wanita 0.120 0.130 0.155
## Tidak Pria 0.065 0.000 0.035
## Wanita 0.120 0.025 0.150
2. Probabilitas Marginal (Marginal Probability)
Probabilitas marginal merupakan penjumlahan dari peluang gabungan untuk masing-masing nilai variabel tertentu, dengan mengabaikan dua variabel lainnya. Hal ini memberikan gambaran menyeluruh terhadap distribusi satu variabel secara independen.
Peluang marginal untuk Sikap = Ya
\[ P(Y = \text{Ya}) = \sum_{Z,X} P(Z,X,Y = \text{Ya}) \]
Sebagai contoh:
\[ P(Y = \text{Ya}) = \frac{11 + 13 + 24 + 24}{200} = \frac{72}{200} = 0.36 \]
Peluang marginal untuk Jenis Kelamin = Pria
\[ P(X = \text{Pria}) = \sum_{Z,Y} P(Z, X = \text{Pria}, Y) \]
Implementasi dalam R
# Hitung probabilitas marginal berdasarkan jenis kelamin (X)
marginal_X <- apply(joint_prob, 2, sum)
marginal_X
## Pria Wanita
## 0.3 0.7
# Hitung probabilitas marginal berdasarkan pengalaman (Z)
marginal_Z <- apply(joint_prob, 1, sum)
marginal_Z
## Ya Tidak
## 0.605 0.395
Tabel Peluang Marginal untuk Jenis Kelamin (X)
| Jenis Kelamin | Probabilitas |
|---|---|
| Pria | 0.545 |
| Wanita | 0.455 |
Tabel Peluang Marginal untuk Pengalaman (Z)
| Pengalaman | Probabilitas |
|---|---|
| Ya | 0.3 |
| Tidak | 0.7 |
Kesimpulan
Analisis probabilitas marginal memberikan informasi tentang persebaran proporsi kategori untuk masing-masing variabel tunggal tanpa mempertimbangkan interaksi dengan variabel lainnya. Ini menjadi pondasi penting sebelum melakukan analisis lanjutan seperti conditional probability atau uji independensi.
3. Probabilitas Bersyarat (Conditional Probability)
Probabilitas bersyarat mengukur peluang suatu kejadian terjadi dengan syarat bahwa kejadian lain telah diketahui. Rumus umumnya:
\[ P(Z \mid X, Y) = \frac{P(Z, X, Y)}{P(X, Y)} \]
Sebagai contoh, misalkan kita ingin mengetahui peluang seseorang pernah mengalami kejadian jika diketahui bahwa dia pria dan memiliki sikap memperdalam ilmu agama:
\[ P(\text{Ya} \mid \text{Pria}, \text{Ya}) = \frac{P(\text{Ya}, \text{Pria}, \text{Ya})}{P(\text{Pria}, \text{Ya})} = \frac{0.055}{0.175} = 0.314 \]
Implementasi dalam R
# Hitung P(X = Pria & Y = Ya)
P_XY <- sum(data_agama_df$P_ZXY[data_agama_df$JenisKelamin == "Pria" & data_agama_df$Sikap == "Ya"])
# Hitung P(Ya | Pria, Ya)
P_Ya_given_PrYa <- data_agama_df$P_ZXY[data_agama_df$Pengalaman == "Ya" &
data_agama_df$JenisKelamin == "Pria" &
data_agama_df$Sikap == "Ya"] / P_XY
P_Ya_given_PrYa
## [1] 0.3142857
##Tabel Peluang Bersyarat
Probabilitas bersyarat dihitung dengan melihat frekuensi pada tiap kombinasi kategori, lalu dibagi dengan jumlah total kategori pengontrol. Ini membantu kita memahami bagaimana peluang sebuah variabel tergantung pada kondisi kategori lain.
Tabel Peluang Bersyarat untuk Pengalaman = Ya
| Jenis Kelamin | Sikap | P(Y |
|---|---|---|
| Pria | Ya | 0.275 |
| Pria | Tidak | 0.300 |
| Pria | Tidak Menentu | 0.425 |
| Wanita | Ya | 0.650 |
| Wanita | Tidak | 0.000 |
| Wanita | Tidak Menentu | 0.350 |
Tabel Peluang Bersyarat untuk Pengalaman = Tidak
| Jenis Kelamin | Sikap | P(Y |
|---|---|---|
| Pria | Ya | 0.235 |
| Pria | Tidak | 0.255 |
| Pria | Tidak Menentu | 0.510 |
| Wanita | Ya | 0.407 |
| Wanita | Tidak | 0.085 |
| Wanita | Tidak Menentu | 0.508 |
Implementasi dalam R
# Buat array 3 dimensi: Pengalaman (Z), Jenis Kelamin (X), Sikap (Y)
data_agama <- array(c(
# Data untuk Sikap = Ya
11, 13, # Pengalaman = Ya: Pria, Wanita
24, 24, # Pengalaman = Tidak: Pria, Wanita
# Data untuk Sikap = Tidak
12, 0,
26, 5,
# Data untuk Sikap = Tidak Menentu
17, 7,
31, 30
),
dim = c(2, 2, 3),
dimnames = list(
Pengalaman = c("Ya", "Tidak"),
JenisKelamin = c("Pria", "Wanita"),
Sikap = c("Ya", "Tidak", "Tidak Menentu")
))
# Hitung peluang bersama
joint_prob <- data_agama / sum(data_agama)
# Hitung peluang bersyarat dengan membagi per baris (setiap Z)
conditional_prob_Ya <- prop.table(data_agama["Ya", , ], margin = 1)
conditional_prob_Tidak <- prop.table(data_agama["Tidak", , ], margin = 1)
# Tampilkan hasil
conditional_prob_Ya
## Sikap
## JenisKelamin Ya Tidak Tidak Menentu
## Pria 0.2750000 0.3000000 0.425000
## Wanita 0.2962963 0.3209877 0.382716
conditional_prob_Tidak
## Sikap
## JenisKelamin Ya Tidak Tidak Menentu
## Pria 0.6500000 0.00000000 0.3500000
## Wanita 0.4067797 0.08474576 0.5084746
Kesimpulan
Pemahaman ini penting untuk menganalisis hubungan antar variabel kategori, khususnya ketika ingin menilai ketergantungan antar variabel berdasarkan konteks yang lebih spesifik.
Conditional independence (kemerdekaan bersyarat) dalam konteks tabel kontingensi merujuk pada kondisi di mana dua variabel menjadi independen satu sama lain setelah mengendalikan variabel ketiga. Dalam istilah probabilitas, variabel X dan Y dikatakan kondisional independen terhadap Z jika:
\[ P(X, Y | Z) = P(X | Z) P(Y | Z) \]
Atau, dalam notasi frekuensi:
\[ \frac{n_{ijk}}{n_{k++}} = \frac{n_{i+k}}{n_{k++}} \times \frac{n_{+jk}}{n_{k++}} \]
Contoh Kasus (Menggunakan Data Agama)
Misalnya kita ingin mengetahui apakah jenis kelamin (X) dan sikap terhadap memperdalam ilmu agama (Y) bersifat independen secara kondisional terhadap pengalaman pernah mengalami kejadian (Z).
Implementasi dalam R
options(kableExtra.latex.load_packages = FALSE)
library(knitr)
## Warning: package 'knitr' was built under R version 4.4.3
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.4.3
# Buat array
data_agama <- array(c(
11, 13, 24, 24, 12, 0, 26, 5, 17, 7, 31, 30
),
dim = c(2, 2, 3),
dimnames = list(
Pengalaman = c("Ya", "Tidak"),
JenisKelamin = c("Pria", "Wanita"),
Sikap = c("Ya", "Tidak", "Tidak Menentu")
))
# Loop dan cetak tabel
for (z in dimnames(data_agama)$Pengalaman) {
df <- as.data.frame.matrix(data_agama[z, , ])
colnames(df) <- c("Ya", "Tidak", "Tidak Menentu")
rownames(df) <- c("Pria", "Wanita")
print(
kable(df, format = "latex", booktabs = TRUE,
caption = paste("Tabel untuk Pengalaman =", z), label = FALSE) %>%
kable_styling(
latex_options = c("hold_position", "striped"),
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)
)
}
Pengujian Conditional Independence
Untuk mengevaluasi apakah sikap dan jenis kelamin saling independen ketika variabel pengalaman dikendalikan, kita bisa melakukan uji chi-square untuk masing-masing strata pengalaman.
Perhitungan Manual untuk Z = Ya
\[ P(\text{Ya} | \text{Pria}, \text{Ya}) = \frac{11}{40} = 0.275 \]
\[ P(\text{Ya} | \text{Wanita}, \text{Ya}) = \frac{13}{20} = 0.65 \]
\[ P(\text{Ya} | \text{Ya}) = \frac{24}{60} = 0.4 \]
Jika nilai \(P(\text{Ya} | \text{Pria}, \text{Ya})\) tidak terlalu berbeda dengan \(P(\text{Ya} | \text{Ya})\), maka bisa dikatakan terdapat conditional independence.
Pengujian dengan R
# Uji chi-square untuk setiap strata pengalaman
chisq_ya <- chisq.test(data_agama["Ya", , ])
chisq_tidak <- chisq.test(data_agama["Tidak", , ])
## Warning in chisq.test(data_agama["Tidak", , ]): Chi-squared approximation may
## be incorrect
chisq_ya
##
## Pearson's Chi-squared test
##
## data: data_agama["Ya", , ]
## X-squared = 0.20023, df = 2, p-value = 0.9047
chisq_tidak
##
## Pearson's Chi-squared test
##
## data: data_agama["Tidak", , ]
## X-squared = 4.3825, df = 2, p-value = 0.1118
Interpretasi Hasil
Uji ini membantu memastikan apakah dua variabel utama saling bergantung satu sama lain setelah memperhitungkan pengaruh variabel ketiga.
Uji Cochran-Mantel-Haenszel (CMH) digunakan untuk menguji apakah dua variabel bersifat independen ketika dikendalikan oleh satu variabel ketiga (variabel pengganggu). Uji ini penting karena:
Prinsip Dasar Uji CMH
Uji ini berangkat dari asumsi bahwa dua variabel utama dapat diuji dalam beberapa strata berdasarkan variabel ketiga. Misalnya, saat menguji pengaruh jenis kelamin dan sikap terhadap agama dengan mempertimbangkan pengalaman (Z), maka dibuat tabel 2x2 untuk masing-masing level Z, lalu dihitung rasio gabungan antar strata.
Jika hubungan antara X dan Y konsisten di semua strata Z, maka diasumsikan ada independensi bersyarat.
Formulasi matematis, hipotesis, dan cara perhitungannya akan dijelaskan dalam bagian berikutnya.
Rumusan Statistik CMH
\[ CMH = \frac{[\sum_k (n_{11k} - \mu_{11k})]^2}{\sum_k \text{var}(n_{11k})} \]
Dengan:
- \(n_{11k}\) = jumlah responden sel
baris 1 kolom 1 pada strata ke-k
- \(\mu_{11k}\) = nilai harapan untuk
sel baris 1 kolom 1 pada strata ke-k, dihitung sebagai:
\[ \mu_{11k} = \frac{n_{1.k} \cdot n_{.1k}}{n_{..k}} \]
\[ \text{var}(n_{11k}) = \frac{n_{1.k} \cdot n_{2.k} \cdot n_{.1k} \cdot n_{.2k}}{n_{..k}^2 (n_{..k} - 1)} \]
Statistik CMH mengikuti distribusi chi-square dengan derajat kebebasan 1.
Tolak \(H_0\) jika nilai CMH > \(\chi^2_{(1)}\) atau nilai p-value < \(\alpha\).
Contoh Penerapan (Data Agama berdasarkan Jenis Kelamin & Sikap, dikontrol oleh Pengalaman)
Misalnya kita ingin menilai apakah hubungan antara jenis kelamin dan sikap memperdalam ilmu agama masih signifikan setelah mengontrol pengalaman. Kita bisa menyusun tabel 2×2 untuk setiap strata Z:
Strata: Pengalaman = Ya
| Ya | Tidak/Tidak Menentu | Total | |
|---|---|---|---|
| Pria | 11 | 29 | 40 |
| Wanita | 13 | 7 | 20 |
Strata: Pengalaman = Tidak
| Ya | Tidak/Tidak Menentu | Total | |
|---|---|---|---|
| Pria | 24 | 57 | 81 |
| Wanita | 24 | 35 | 59 |
Perhitungan Manual CMH
1. Menghitung nilai harapan \(\mu_{11k}\)
Rumus nilai harapan \(\mu_{11k}\) untuk setiap strata \(k\):
\[ \mu_{11k} = \frac{n_{1{+}k} \times n_{+1k}}{n_{++k}} \]
Keterangan:
- \(n_{1{+}k}\): jumlah Pria dalam
strata pengalaman \(k\)
- \(n_{+1k}\): jumlah “Ya” memperdalam
ilmu agama dalam strata pengalaman \(k\)
- \(n_{++k}\): total individu dalam
strata \(k\)
Untuk strata “Pernah Mengalami” (\(k = 1\)):
\[ \mu_{111} = \frac{40 \times 24}{60} = \frac{960}{60} = 16.00 \]
Untuk strata “Tidak Pernah Mengalami” (\(k = 2\)):
\[ \mu_{112} = \frac{81 \times 48}{140} = \frac{3888}{140} \approx 27.77 \]
2. Menghitung Varians \(\text{Var}(n_{11k})\)
Rumus:
\[ \text{Var}(n_{11k}) = \frac{n_{1{+}k} n_{0{+}k} n_{+1k} n_{+0k}}{n_{++k}^2 (n_{++k} - 1)} \]
Untuk \(k = 1\):
\[ \text{Var}(n_{111}) = \frac{40 \times 20 \times 24 \times 36}{60^2 (60 - 1)} = \frac{691200}{212400} \approx 9.36 \]
Untuk \(k = 2\):
\[ \text{Var}(n_{112}) = \frac{81 \times 59 \times 48 \times 92}{140^2 (140 - 1)} = \frac{21299088}{2707600} \approx 7.87 \]
3. Menghitung Statistik CMH
\[ X^2_{CMH} = \frac{\left(\sum_k (n_{11k} - \mu_{11k})\right)^2}{\sum_k \text{Var}(n_{11k})} \]
\[ X^2_{CMH} = \frac{(11 - 16 + 24 - 27.77)^2}{9.36 + 7.87} = \frac{(-8.77)^2}{17.23} = \frac{76.94}{17.23} \approx 4.47 \]
Bandingkan dengan nilai kritis \(\chi^2_{0.05,1} = 3.841\)
Karena \(4.47 > 3.841\), maka hipotesis nol ditolak. Artinya, ada hubungan signifikan antara jenis kelamin dan sikap memperdalam ilmu agama setelah dikontrol oleh pengalaman.
Implementasi Uji CMH di R
library(vcdExtra)
## Warning: package 'vcdExtra' was built under R version 4.4.3
## Loading required package: vcd
## Warning: package 'vcd' was built under R version 4.4.3
## Loading required package: grid
## Loading required package: gnm
## Warning: package 'gnm' was built under R version 4.4.3
# Bentuk array CMH berdasarkan data sebelumnya
data_cmh <- array(c(11, 29, 13, 7, # strata Ya (Pernah Mengalami)
24, 57, 24, 35), # strata Tidak Pernah
dim = c(2, 2, 2),
dimnames = list(
JenisKelamin = c("Pria", "Wanita"),
Sikap = c("Ya", "Tidak/Tidak Menentu"),
Pengalaman = c("Ya", "Tidak")
))
# Uji CMH
mantelhaen.test(data_cmh, correct = FALSE)
##
## Mantel-Haenszel chi-squared test without continuity correction
##
## data: data_cmh
## Mantel-Haenszel X-squared = 6.994, df = 1, p-value = 0.008178
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.2509041 0.8202474
## sample estimates:
## common odds ratio
## 0.4536556
Kesimpulan
Jika nilai statistik CMH lebih besar dari nilai kritis \(\chi^2\) atau jika p-value < 0.05, maka kita menolak hipotesis nol. Ini berarti ada hubungan signifikan antara jenis kelamin dan sikap memperdalam ilmu agama setelah dikontrol oleh pengalaman.
Hasil ini menunjukkan pentingnya mempertimbangkan variabel kontrol dalam analisis tabel kontingensi tiga arah.
Dalam analisis data kontingensi 2×2×K, kita sering menemukan sejumlah tabel parsial berdasarkan stratifikasi variabel kontrol. Bila nilai-nilai odds ratio kondisional dari tiap strata cenderung seragam, kita dapat menghitung satu nilai ringkasan yang disebut sebagai odds ratio bersama.
Penaksiran Odds Ratio Bersama
Odds ratio bersama dapat dihitung menggunakan estimasi Mantel–Haenszel dengan rumus:
\[ \hat{\theta}_{MH} = \frac{\sum_{k=1}^K \left( \frac{n_{11k}n_{22k}}{n_{++k}} \right)}{\sum_{k=1}^K \left( \frac{n_{12k}n_{21k}}{n_{++k}} \right)} \]
Dengan:
- \(n_{11k}\): jumlah Pria yang
menjawab “Ya” dalam strata pengalaman \(k\)
- \(n_{12k}\): Pria yang menjawab
“Tidak/Tidak Menentu” - \(n_{21k}\):
Wanita yang menjawab “Ya”
- \(n_{22k}\): Wanita yang menjawab
“Tidak/Tidak Menentu” - \(n_{++k}\):
total responden dalam strata \(k\)
Dari data kita:
\[ \hat{\theta}_{MH} = \frac{(11 \times 7)/60 + (24 \times 35)/140}{(29 \times 13)/60 + (57 \times 24)/140} = \frac{(77/60) + (840/140)}{(377/60) + (1368/140)} = \frac{1.283 + 6.00}{6.283 + 9.77} = \frac{7.283}{16.05} = 0.454 \]
Standard Error log Odds Ratio Bersama
Standard error log odds ratio bersama dapat dihitung menggunakan rumus:
\[ \hat{\sigma}^2[\log(\hat{\theta}_{MH})] = \frac{\sum (n_{11k} + n_{12k})(n_{11k} + n_{21k})}{2 (\sum n_{11k}n_{22k}/n_{++k})^2} + \frac{\sum (n_{21k} + n_{22k})(n_{12k} + n_{22k})}{2 (\sum n_{12k}n_{21k}/n_{++k})^2} \]
Untuk keperluan praktik, hasil dari fungsi
mantelhaen.test() dalam R akan langsung memberikan nilai
log odds ratio, standard error, dan interval kepercayaan.
Contoh Odds Ratio Bersama
Berikut contoh hasil odds ratio bersama menggunakan data memperdalam ilmu agama:
\[ \hat{\theta}_{MH} = \frac{(11 \times 7)/60 + (24 \times 35)/140}{(29 \times 13)/60 + (57 \times 24)/140} = \frac{7.283}{16.05} = 0.454 \]
\[ SE(\log \hat{\theta}_{MH}) = 0.25 \]
\[ \log(0.454) \pm 1.96 \times 0.25 = -0.791 \pm 0.49 \]
\[ 1.98 \leq \theta \leq 2.38 \]
\[ \text{CI OR} = (e^{-1.281}, e^{-0.301}) = (0.28, 0.74) \]
Interpretasi
Kesimpulan
Odds ratio bersama digunakan untuk merangkum hubungan antara dua variabel kategori utama dalam kehadiran variabel kontrol. Pendekatan Mantel–Haenszel memungkinkan perhitungan odds ratio yang memperhitungkan stratifikasi. Dalam kasus ini, kita menemukan bahwa jenis kelamin masih berhubungan secara signifikan dengan sikap memperdalam ilmu agama setelah dikontrol oleh pengalaman pribadi.
Implementasi di R
Kita dapat menghitung odds ratio bersama, standard error, dan
interval kepercayaan menggunakan paket epitools dalam
R:
# Paket yang dibutuhkan
library(epitools)
##
## Attaching package: 'epitools'
## The following object is masked from 'package:vcdExtra':
##
## expand.table
## The following object is masked from 'package:vcd':
##
## oddsratio
# Susun array 2x2x2 berdasarkan pengalaman
agama_array <- array(c(11, 29, 13, 7, # strata Pengalaman = Ya
24, 57, 24, 35), # strata Pengalaman = Tidak
dim = c(2, 2, 2),
dimnames = list(
JenisKelamin = c("Pria", "Wanita"),
Sikap = c("Ya", "Tidak/Tidak Menentu"),
Pengalaman = c("Ya", "Tidak")
))
# Hitung odds ratio bersama
mh_result <- mantelhaen.test(agama_array, correct = FALSE)
print(mh_result)
##
## Mantel-Haenszel chi-squared test without continuity correction
##
## data: agama_array
## Mantel-Haenszel X-squared = 6.994, df = 1, p-value = 0.008178
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.2509041 0.8202474
## sample estimates:
## common odds ratio
## 0.4536556
Output dari fungsi ini akan mencakup estimasi odds ratio bersama, nilai statistik uji, p-value, serta confidence interval.
Uji Homogenitas Odds Ratio dengan Statistik Breslow-Day
Homogenitas Asosiasi dalam Tabel Kontingensi Tiga Arah
Definisi Asosiasi Homogen
Disebut asosiasi homogen jika nilai odds ratio untuk tiap tabel parsial bernilai sama:
\(\theta_{xy(1)} = \theta_{xy(2)} = \cdots = \theta_{xy(k)}\)
Bila odds ratio tersebut konsisten di semua strata untuk variabel kontrol \(Z\), artinya tidak ada interaksi antara variabel \(X\) dan \(Y\) dalam memengaruhi \(Z\).
Namun, jika odds ratio bervariasi di antara level \(Z\), maka ada interaksi antara \(X\) dan \(Y\) dalam hubungannya terhadap \(Z\).
Pengujian Homogenitas dengan Statistik Breslow-Day
Hipotesis:
Statistik Uji Breslow-Day (BD) digunakan untuk menguji apakah odds ratio antar strata tersebut homogen:
\[ X^2_{HBD} = \sum_{j=1}^K \frac{(a_j - \tilde{a}_j)^2}{\text{Var}(a_j|\hat{OR}_{MH})} \]
dengan:
- \(a_j\) adalah jumlah kasus terpapar
yang diamati dalam strata ke-\(j\)
- \(\tilde{a}_j\) adalah nilai yang
diharapkan (ekspektasi) dari kasus terpapar
- \(\text{Var}(a_j|\hat{OR}_{MH})\)
adalah varians dari nilai ekspektasi tersebut
Langkah-Langkah Perhitungan
Menghitung Rasio Odds Gabungan (Mantel–Haenszel OR)
Estimasi odds ratio gabungan dihitung sebagai:
\[ \hat{OR}_{MH} = \frac{\sum_{j=1}^K \frac{a_j d_j}{n_j}}{\sum_{j=1}^K \frac{b_j c_j}{n_j}} \]
dengan \(a_j, b_j, c_j, d_j\) adalah elemen pada tabel 2x2 di strata \(j\) dan \(n_j\) total observasi di strata \(j\)
Menentukan Nilai Ekspektasi \(\tilde{a}_j\)
Nilai ini diperoleh dengan menyelesaikan persamaan kuadrat:
\[ -m_{1j}n_{1j}\hat{OR}_{MH} + ((n_{2j} - m_{1j}) + \hat{OR}_{MH}(n_{1j} + m_{1j}))\tilde{a}_j + (1 - \hat{OR}_{MH})\tilde{a}_j^2 = 0 \]
dengan \(m_{1j}\) adalah total terpapar di strata \(j\), \(n_{1j}\) adalah total kasus di strata \(j\) Solusi kuadrat dengan syarat \(0 < \tilde{a}_j \leq \min(n_{1j}, m_{1j})\)
Menghitung Varians \(\text{Var}(a_j|\hat{OR}_{MH})\)
Varians dihitung sebagai:
\[ \text{Var}(a_j|\hat{OR}_{MH}) = \left(\frac{1}{\tilde{a}_j} + \frac{1}{\tilde{b}_j} + \frac{1}{\tilde{c}_j} + \frac{1}{\tilde{d}_j}\right)^{-1} \]
dengan:
Menghitung Statistik Breslow-Day
Gunakan rumus:
\[ X^2_{HBD} = \sum_{j=1}^K \frac{(a_j - \tilde{a}_j)^2}{\text{Var}(a_j|\hat{OR}_{MH})} \]
Koreksi Tarone
Koreksi ini digunakan untuk memperbaiki bias dengan rumus:
\[ X^2_{HBDT} = X^2_{HBD} - \frac{\left(\sum_{j=1}^K (a_j - \tilde{a}_j)\right)^2}{\sum_{j=1}^K \text{Var}(a_j|\hat{OR}_{MH})} \]
P-value dan Keputusan Hipotesis
Kesimpulan - Statistik Breslow-Day
Uji \(X^2_{HBD}\) bertujuan mengetahui apakah rasio odds di tiap strata (misalnya untuk data memperdalam ilmu agama) homogen atau tidak. Nilai ekspektasi \(\tilde{a}_j\) diperoleh melalui penyelesaian kuadrat. Koreksi Tarone \(X^2_{HBDT}\) digunakan agar hasilnya tidak bias. Terakhir, p-value dihitung menggunakan distribusi chi-square untuk memutuskan apakah hipotesis homogenitas (\(H_0\)) dapat ditolak atau tidak.
Ilustrasi Perhitungan Manual Menggunakan Data Memperdalam Ilmu Agama
Data disusun menjadi dua strata berdasarkan variabel Z (Pernah Mengalami Kejadian):
| Z (Pengalaman) | Jenis Kelamin | Y = Ya | Y = Tidak/Tidak Menentu | Total |
|---|---|---|---|---|
| Ya | Pria | 11 | 29 | 40 |
| Ya | Wanita | 13 | 7 | 20 |
| Tidak | Pria | 24 | 57 | 81 |
| Tidak | Wanita | 24 | 35 | 59 |
Langkah-Langkah:
Kesimpulan: Karena \(p > 0.05\), maka tidak ada bukti cukup untuk menolak hipotesis bahwa odds ratio memperdalam ilmu agama homogen di seluruh strata.
# Fungsi untuk Uji Breslow-Day dengan koreksi Tarone
breslowday.test <- function(x) {
# Estimasi common OR menggunakan Mantel-Haenszel
or.hat.mh <- mantelhaen.test(x)$estimate
# Banyaknya strata
K <- dim(x)[3]
# Inisialisasi nilai statistik
X2.HBD <- 0
a <- tildea <- Var.a <- numeric(K)
for (j in 1:K) {
# Marginal baris dan kolom pada strata ke-j
mj <- apply(x[,,j], 1, sum) # baris
nj <- apply(x[,,j], 2, sum) # kolom
# Solusi tilde_a_j (expected frek)
coef <- c(-mj[1]*nj[1] * or.hat.mh, nj[2] - mj[1] + or.hat.mh*(nj[1] + mj[1]), 1 - or.hat.mh)
sols <- Re(polyroot(coef))
tildeaj <- sols[(0 < sols) & (sols <= min(nj[1], mj[1]))]
aj <- x[1,1,j]
# Sel ekspektasi lainnya
tildebj <- mj[1] - tildeaj
tildecj <- nj[1] - tildeaj
tildedj <- mj[2] - tildecj
# Varian a_j
Var.aj <- (1/tildeaj + 1/tildebj + 1/tildecj + 1/tildedj)^(-1)
# Kontribusi ke statistik Breslow-Day
X2.HBD <- X2.HBD + as.numeric((aj - tildeaj)^2 / Var.aj)
a[j] <- aj
tildea[j] <- tildeaj
Var.a[j] <- Var.aj
}
# Koreksi Tarone
X2.HBDT <- as.numeric(X2.HBD - (sum(a) - sum(tildea))^2 / sum(Var.a))
p <- 1 - pchisq(X2.HBDT, df = K - 1)
res <- list(X2.HBD = X2.HBD, X2.HBDT = X2.HBDT, p = p)
class(res) <- "bdtest"
return(res)
}
# Fungsi print hasil uji
print.bdtest <- function(x) {
cat("Breslow-Day Test (Tarone correction):\n")
cat("Breslow-Day Chi-squared =", x$X2.HBD, "\n")
cat("Tarone-corrected Chi-squared =", x$X2.HBDT, "\n")
cat("P-value =", x$p, "\n")
}
breslowday.test(agama_array)
## Breslow-Day Test (Tarone correction):
## Breslow-Day Chi-squared = 2.584873
## Tarone-corrected Chi-squared = 2.58338
## P-value = 0.1079908
# Instal dan muat paket jika belum ada
if (!require(DescTools)) install.packages("DescTools")
## Loading required package: DescTools
## Warning: package 'DescTools' was built under R version 4.4.3
library(DescTools)
# Array data 2x2x2 berdasarkan pengalaman
agama_array <- array(c(11, 29, 13, 7, # strata Pengalaman = Ya
24, 57, 24, 35), # strata Pengalaman = Tidak
dim = c(2, 2, 2),
dimnames = list(
JenisKelamin = c("Pria", "Wanita"),
Sikap = c("Ya", "Tidak/Tidak Menentu"),
Pengalaman = c("Ya", "Tidak")
))
# Uji Breslow-Day homogenitas OR antar strata
breslow_test <- BreslowDayTest(agama_array)
print(breslow_test)
##
## Breslow-Day test on Homogeneity of Odds Ratios
##
## data: agama_array
## X-squared = 2.5849, df = 1, p-value = 0.1079
Interpretasi Hasil
Generalized Linear Model (GLM) merupakan pengembangan dari model regresi linier klasik yang memperluas kemampuannya dalam menangani variabel respons yang tidak mengikuti distribusi normal serta hubungan non-linear antara variabel prediktor dan ekspektasi respons.
GLM terdiri dari tiga komponen utama:
Sebuah distribusi termasuk ke dalam keluarga exponential family apabila memiliki bentuk fungsi probabilitas:
\[ f(y; \theta, \phi) = \exp\left\{ \frac{y \theta - b(\theta)}{a(\phi)} + c(y, \phi) \right\} \]
Keterangan:
Distribusi yang termasuk dalam keluarga exponential family meliputi:
Contoh Transformasi: Distribusi Binomial
Distribusi Binomial dengan fungsi:
\[ P(Y = y) = \binom{n}{y} \pi^y (1 - \pi)^{n - y} \]
Dapat dituliskan dalam bentuk exponential family sebagai:
\[ P(Y = y) = \exp \left\{ \log\binom{n}{y} + y \log\left( \frac{\pi}{1 - \pi} \right) + n \log(1 - \pi) \right\} \]
Dengan substitusi:
Sehingga fungsi binomial dapat dituliskan dalam format exponential family.
Model regresi logistik digunakan ketika variabel respons bersifat biner (hanya memiliki dua kategori: 0 dan 1). Model ini memanfaatkan fungsi aktivasi sigmoid untuk membatasi nilai output dalam rentang 0 hingga 1, membentuk kurva berbentuk huruf S (S-shaped).
Regresi logistik memodelkan probabilitas keanggotaan suatu observasi dalam satu dari dua kategori. Model ini sangat berguna untuk klasifikasi, seperti memprediksi apakah seseorang memiliki penyakit tertentu, atau apakah email tergolong spam.
Contoh aplikasi umum:
Keunggulan Regresi Logistik
Jika output dari fungsi sigmoid:
Dengan demikian, regresi logistik sangat relevan untuk prediksi kategori biner berbasis data input numerik atau kategorik.
Fungsi Sigmoid
Fungsi sigmoid didefinisikan sebagai:
\[ f(x) = \frac{1}{1 + e^{-x}} \]
Jika hasil prediksi lebih besar dari ambang batas (misalnya 0.5), maka hasil diklasifikasikan sebagai 1. Jika di bawah 0.5, maka diklasifikasikan sebagai 0.
Simulasi dan Visualisasi Regresi Logistik
# Simulasi data untuk regresi logistik
set.seed(42)
n <- 100
x <- seq(-4, 4, length.out = n)
log_odds <- -0.5 + 1.5 * x
prob <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, prob)
# Buat data frame
data <- data.frame(x = x, y = y, prob = prob)
# Visualisasi menggunakan base R
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 (0/1)", "Kurva Logistik", "Ambang 0.5"),
col = c("gray60", "blue", "red"),
pch = c(16, NA, NA),
lty = c(NA, 1, 2),
lwd = c(NA, 2, 1),
pt.cex = 1.5,
bty = "n")
Plot Probabilitas Prediksi
Berdasarkan hasil prediksi model, kita dapat memvisualisasikan probabilitas sebagai berikut:
model <- glm(y ~ x, data = data, family = binomial)
data$pred <- predict(model, type = "response")
plot(data$x, data$y, pch = 16, col = "grey")
curve(predict(model, newdata = data.frame(x = x), type = "response")[order(x)],
add = TRUE, col = "blue", lwd = 2)
Spesifikasi Model Regresi Logistik
Kurva sigmoid dalam regresi logistik memperlihatkan hubungan non-linier antara prediktor dan probabilitas hasil.
Fungsi link (logit): \[ g(\mu) = \log\left( \frac{\mu}{1 - \mu} \right) \]
Model: \[ \log\left( \frac{\mu}{1 - \mu} \right) = X \beta \]
Fungsi inverse link: \[ \mu = \frac{\exp(X\beta)}{1 + \exp(X\beta)} \]
Regresi logistik pada GLM umumnya menggunakan metode Maximum Likelihood Estimation (MLE). Fungsi log-likelihood:
\[ l(\beta) = \sum_{i=1}^n \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right] \]
Dengan:
\[ \pi_i = \frac{\exp(x_i^T \beta)}{1 + \exp(x_i^T \beta)} \]
Estimasi parameter diperoleh melalui metode numerik seperti Newton-Raphson atau Fisher Scoring.
Contoh Estimasi di R
set.seed(123)
n <- 200
x <- rnorm(n)
p <- 1 / (1 + exp(-(-0.5 + 2*x)))
y <- rbinom(n, 1, p)
data <- data.frame(y, x)
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.8374 0.1954 -4.286 1.82e-05 ***
## x 2.0262 0.3016 6.718 1.84e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 261.37 on 199 degrees of freedom
## Residual deviance: 177.26 on 198 degrees of freedom
## AIC: 181.26
##
## Number of Fisher Scoring iterations: 5
Estimasi ini memberikan koefisien regresi dan ukuran signifikansi untuk memeriksa pengaruh prediktor terhadap probabilitas respons.
GLM adalah kerangka kerja yang fleksibel untuk menganalisis berbagai
jenis data, baik numerik maupun kategorik. Regresi logistik adalah salah
satu contoh paling umum dari GLM, sangat bermanfaat dalam klasifikasi
biner. Estimasi parameter dilakukan dengan MLE dan dapat
diimplementasikan secara efisien dengan fungsi glm di
R.
Contoh 3
Contoh ini menggunakan dataset iris dari R base package,
namun dikonversi menjadi kasus klasifikasi biner. Kita ingin memodelkan
apakah bunga termasuk spesies setosa (1) atau bukan (0),
berdasarkan variabel Sepal.Length dan
Petal.Length.
Persiapan Data
data(iris)
iris$target <- ifelse(iris$Species == "setosa", 1, 0)
summary(iris[, c("target", "Sepal.Length", "Petal.Length")])
## target Sepal.Length Petal.Length
## Min. :0.0000 Min. :4.300 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:5.100 1st Qu.:1.600
## Median :0.0000 Median :5.800 Median :4.350
## Mean :0.3333 Mean :5.843 Mean :3.758
## 3rd Qu.:1.0000 3rd Qu.:6.400 3rd Qu.:5.100
## Max. :1.0000 Max. :7.900 Max. :6.900
Estimasi Model Regresi Logistik
model_real <- glm(target ~ Sepal.Length + Petal.Length, data = iris, family = binomial)
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(model_real)
##
## Call:
## glm(formula = target ~ Sepal.Length + Petal.Length, family = binomial,
## data = iris)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -78.26 316177.72 0.000 1.000
## Sepal.Length 39.83 81738.89 0.000 1.000
## Petal.Length -48.60 43659.68 -0.001 0.999
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1.9095e+02 on 149 degrees of freedom
## Residual deviance: 5.0543e-09 on 147 degrees of freedom
## AIC: 6
##
## Number of Fisher Scoring iterations: 25
Interpretasi Koefisien
setosa saat panjang
sepal dan petal bernilai nolsetosasetosaOdds Ratio:
exp(coef(model_real))
## (Intercept) Sepal.Length Petal.Length
## 1.023812e-34 1.993132e+17 7.814816e-22
setosasetosaVisualisasi Prediksi
iris$prob <- predict(model_real, type = "response")
iris$pred_class <- ifelse(iris$prob > 0.5, "setosa", "non-setosa")
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
ggplot(iris, aes(x = Sepal.Length, y = Petal.Length, color = pred_class, shape = Species)) +
geom_point(size = 3) +
labs(title = "Prediksi Regresi Logistik: Klasifikasi Setosa",
x = "Panjang Sepal (cm)",
y = "Panjang Petal (cm)",
color = "Prediksi", shape = "Fakta") +
theme_minimal()
Evaluasi Model
# Confusion matrix
table(Predicted = iris$pred_class, Actual = iris$Species == "setosa")
## Actual
## Predicted FALSE TRUE
## non-setosa 100 0
## setosa 0 50
Interpretasi:
Kesimpulan
setosa dari panjang
sepal dan petalsetosasetosaModel regresi Poisson digunakan untuk memodelkan data cacah (count data), yaitu data yang bernilai bulat non-negatif. Model ini merupakan bagian dari GLM yang mengasumsikan bahwa distribusi dari variabel respons mengikuti distribusi Poisson.
Distribusi Poisson memiliki fungsi probabilitas:
\[ P(Y = y) = \frac{e^{-\lambda} \lambda^y}{y!} \]
Fungsi tersebut dapat dituliskan dalam format exponential family:
\[ f(y; \theta) = \exp \left\{ y \log(\lambda) - \lambda - \log(y!) \right\} \]
Dengan substitusi:
Sehingga distribusi Poisson termasuk dalam exponential family.
Fungsi link kanonik untuk Poisson adalah fungsi logaritma:
\[ g(\mu) = \log(\mu) \]
Model linier menjadi:
\[ \log(\mu_i) = x_i^T \beta \]
Fungsi inverse link:
\[ \mu_i = \exp(x_i^T \beta) \]
Estimasi parameter \(\beta\) dilakukan menggunakan Maximum Likelihood Estimation (MLE). Fungsi log-likelihood:
\[ l(\beta) = \sum_{i=1}^n \left[ y_i x_i^T \beta - \exp(x_i^T \beta) - \log(y_i!) \right] \]
Parameter dapat diestimasi menggunakan algoritma numerik seperti Newton-Raphson.
Contoh 1
Simulasi data untuk regresi Poisson:
set.seed(42)
n <- 200
x <- rnorm(n)
lambda <- exp(0.3 + 0.6 * x)
y <- rpois(n, lambda)
data <- data.frame(y, x)
# Estimasi model
poisson_model <- glm(y ~ x, data = data, family = poisson)
summary(poisson_model)
##
## Call:
## glm(formula = y ~ x, family = poisson, data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.21817 0.06712 3.250 0.00115 **
## x 0.58748 0.06288 9.343 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 333.72 on 199 degrees of freedom
## Residual deviance: 244.46 on 198 degrees of freedom
## AIC: 587.36
##
## Number of Fisher Scoring iterations: 5
Model ini menghasilkan koefisien regresi yang menjelaskan pengaruh variabel prediktor terhadap jumlah kejadian rata-rata (mean count).
Plot Hasil Prediksi
Visualisasi prediksi hasil model regresi Poisson dapat dilakukan seperti berikut:
plot(data$x, data$y, pch = 16, col = "darkgray", main = "Data dan Hasil Prediksi")
newdata <- data.frame(x = seq(min(x), max(x), length.out = 100))
pred <- predict(poisson_model, newdata = newdata, type = "response")
lines(newdata$x, pred, col = "blue", lwd = 2)
Diagnostik dan Overdispersion
Salah satu asumsi dari model Poisson adalah kesamaan antara nilai ekspektasi dan varians (\(E[Y] = Var[Y]\)). Bila varians melebihi nilai rata-ratanya, maka terjadi kondisi overdispersion.
Deteksi overdispersion dapat dilakukan dengan:
dispersion <- sum(residuals(poisson_model, type = "pearson")^2) / poisson_model$df.residual
dispersion
## [1] 1.144237
Jika nilai hasilnya > 1, maka overdispersion mungkin terjadi dan kita bisa mempertimbangkan model alternatif seperti regresi binomial negatif.
Model Poisson sangat bermanfaat dalam analisis data cacah, tetapi penting untuk selalu mengevaluasi potensi overdispersion.
Contoh 2
Sebagai contoh lain, kita menggunakan dataset
InsectSprays dari R base package. Dataset ini mencatat
jumlah serangga yang masih hidup setelah diberikan perlakuan dengan enam
jenis insektisida berbeda.
Deskripsi Dataset:
count: jumlah serangga yang tersisa (variabel
target)spray: jenis insektisida (A hingga F)data("InsectSprays")
head(InsectSprays)
summary(InsectSprays)
## count spray
## Min. : 0.00 A:12
## 1st Qu.: 3.00 B:12
## Median : 7.00 C:12
## Mean : 9.50 D:12
## 3rd Qu.:14.25 E:12
## Max. :26.00 F:12
Pemodelan Regresi Poisson
poisson_model <- glm(count ~ spray, data = InsectSprays, family = poisson)
summary(poisson_model)
##
## Call:
## glm(formula = count ~ spray, family = poisson, data = InsectSprays)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.67415 0.07581 35.274 < 2e-16 ***
## sprayB 0.05588 0.10574 0.528 0.597
## sprayC -1.94018 0.21389 -9.071 < 2e-16 ***
## sprayD -1.08152 0.15065 -7.179 7.03e-13 ***
## sprayE -1.42139 0.17192 -8.268 < 2e-16 ***
## sprayF 0.13926 0.10367 1.343 0.179
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 409.041 on 71 degrees of freedom
## Residual deviance: 98.329 on 66 degrees of freedom
## AIC: 376.59
##
## Number of Fisher Scoring iterations: 5
Interpretasi Koefisien
exp(coef(poisson_model))
## (Intercept) sprayB sprayC sprayD sprayE sprayF
## 14.5000000 1.0574713 0.1436782 0.3390805 0.2413793 1.1494253
Visualisasi Prediksi
InsectSprays$predicted <- predict(poisson_model, type = "response")
library(ggplot2)
ggplot(InsectSprays, aes(x = spray, y = count, color = spray)) +
geom_jitter(width = 0.2, alpha = 0.6) +
geom_point(aes(y = predicted), shape = 18, size = 3, color = "black") +
labs(title = "Prediksi Jumlah Serangga berdasarkan Jenis Insektisida",
x = "Jenis Spray",
y = "Jumlah Serangga") +
theme_minimal()
Evaluasi Model
# Plot residual
plot(poisson_model$residuals, main = "Residual Plot", ylab = "Residual", xlab = "Index", pch = 19, col = "blue")
abline(h = 0, col = "red", lty = 2)
Kesimpulan Model Poisson
Dalam Generalized Linear Model (GLM), proses inferensial membutuhkan pemahaman mendalam terhadap ekspektasi dan varians dari estimasi parameter. Hal ini sangat penting dalam membangun uji statistik seperti uji Wald, Likelihood Ratio test, dan interval kepercayaan.
Sub-bab ini mengulas dua konsep penting dalam inferensi GLM, yaitu ekspektasi dan varians dari estimator parameter. Pemahaman terhadap ekspektasi dan varians sangat penting untuk menilai ketepatan dan keandalan hasil estimasi.
1. Ekspektasi Estimator
Ekspektasi digunakan untuk menilai apakah suatu estimator bersifat tidak bias. Sebuah estimator \(\hat{\beta}\) dikatakan unbiased jika:
\[ \mathbb{E}[\hat{\beta}] = \beta \]
Dalam konteks GLM, estimator Maximum Likelihood (MLE) dari \(\hat{\beta}\) memiliki sifat tidak bias secara asimtotik (asymptotically unbiased).
2. Varians Estimator
Varians mengukur tingkat presisi dari estimator. Untuk GLM, varians dari \(\hat{\beta}\) dapat dihampiri dengan:
\[ \text{Var}(\hat{\beta}) \approx (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \]
di mana \(\mathbf{W}\) adalah matriks bobot yang ditentukan oleh bentuk distribusi dan fungsi link yang digunakan.
3. Distribusi Asimptotik Estimator
Bagian ini membahas distribusi dari estimator parameter ketika ukuran sampel cukup besar. Distribusi asimptotik digunakan sebagai dasar untuk uji statistik dan penarikan kesimpulan.
Ketika ukuran sampel besar, distribusi dari \(\hat{\beta}\) mendekati distribusi normal:
\[ \hat{\beta} \sim \mathcal{N}(\beta, \text{Var}(\hat{\beta})) \]
Distribusi ini digunakan dalam penyusunan: - Uji Wald - Interval kepercayaan - Perhitungan p-value
4. Varians Tidak Konstan dalam GLM
Sub-bab ini menjelaskan bahwa varians dalam GLM tidak konstan, tidak seperti regresi linear klasik. Varians dalam GLM bergantung pada nilai ekspektasi dari variabel respons dan bentuk distribusi.
Tidak seperti model OLS yang mengasumsikan varian konstan (homoskedastisitas):
\[ \text{Var}(Y_i) = \sigma^2 \]
Dalam GLM:
\[ \text{Var}(Y_i) = \phi V(\mu_i) \]
di mana:
- \(\phi\): parameter dispersi
- \(V(\mu)\): fungsi varians spesifik
untuk setiap distribusi
Contoh fungsi varians:
- Distribusi Poisson: \(V(\mu) =
\mu\)
- Distribusi Binomial: \(V(\mu) = \mu(1 -
\mu)\)
Simulasi Regresi Poisson di R
Simulasi ini bertujuan untuk memberikan ilustrasi praktis penerapan regresi Poisson dalam GLM menggunakan bahasa pemrograman R.
# Simulasi data
set.seed(2025)
x <- rnorm(100)
mu <- exp(0.7 + 0.9 * x)
y <- rpois(100, mu)
model <- glm(y ~ x, family = poisson)
summary(model)
##
## Call:
## glm(formula = y ~ x, family = poisson)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.63277 0.08031 7.879 3.29e-15 ***
## x 0.94733 0.05585 16.964 < 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: 391.39 on 99 degrees of freedom
## Residual deviance: 102.76 on 98 degrees of freedom
## AIC: 334.66
##
## Number of Fisher Scoring iterations: 5
Kesimpulan Penting:
Bagian ini merangkum poin-poin utama yang dibahas mengenai ekspektasi dan varians dalam GLM.
Sub-bab ini menjelaskan bagaimana cara menurunkan ekspektasi dan varians dalam GLM secara matematis, khususnya melalui pendekatan fungsi momen dan formulasi distribusi eksponensial.
Ekspektasi
Ekspektasi dapat diturunkan menggunakan fungsi momen:
\[ \mathbb{E}(Y) = \int y f(y; \theta) dy = \mu \]
Untuk distribusi dari keluarga eksponensial:
\[ \log f(y; \theta) = a(y) + b(\theta)y + c(\theta) \quad \text{atau} \quad \log f(y; \theta) = y\theta - b(\theta) + c(y) \]
Turunan pertama dari log-likelihood menghasilkan fungsi skor:
\[ U(\theta) = \frac{\partial \ell}{\partial \theta} = y - b'(\theta) \]
Dan nilai ekspektasi dari turunan pertama ini:
\[ \mathbb{E}[U(\theta)] = \mathbb{E}[y - b'(\theta)] = \mu - b'(\theta) = 0 \]
Maka:
\[ \mu = b'(\theta) \]
Varians
Sub-bab ini melanjutkan dari penurunan ekspektasi menuju ke turunan kedua fungsi log-likelihood, yang digunakan untuk menurunkan bentuk varian.
Turunan kedua dari log-likelihood:
\[ \frac{\partial^2 \ell}{\partial \theta^2} = -b''(\theta) \]
Dengan demikian:
\[ \text{Var}(Y) = b''(\theta) = \phi V(\mu) \]
Bab ini membahas pendekatan yang digunakan untuk mengestimasi parameter dalam GLM, dengan fokus pada metode Maximum Likelihood dan teknik numerik seperti Newton-Raphson.
1. Maximum Likelihood Estimation (MLE)
MLE bertujuan untuk menentukan nilai parameter yang memaksimalkan fungsi likelihood atau log-likelihood. Syarat optimum:
Namun, karena bentuk GLM sering tidak memiliki solusi eksplisit, diperlukan metode numerik.
2. Metode Optimisasi Newton-Raphson
Newton-Raphson menggunakan pendekatan numerik berbasis turunan:
- Score vector \(U(\beta)\)
- Hessian matrix \(H(\beta)\)
Iterasi:
\[ \beta^{(t+1)} = \beta^{(t)} - H^{-1}(\beta^{(t)}) U(\beta^{(t)}) \]
3. Fisher Scoring
Metode ini adalah modifikasi dari Newton-Raphson, di mana matriks Hessian diganti dengan matriks informasi Fisher. Hal ini sering lebih stabil secara numerik.
4. IRLS (Iteratively Reweighted Least Square)
Merupakan pendekatan alternatif yang juga berbasis iterasi. IRLS sangat umum digunakan dalam implementasi praktis GLM karena efisien dan mudah diadaptasi.
5. Implementasi Newton-Raphson
Untuk setiap parameter \(\beta_j\), score function dan turunan kedua dirumuskan sebagai:
\[ U_j(\beta) = \frac{\partial \log L(\beta)}{\partial \beta_j} \] \[ H_{jk}(\beta) = \frac{\partial^2 \log L(\beta)}{\partial \beta_j \partial \beta_k} \]
Dengan pendekatan Taylor:
\[ U(\beta^*) \approx U(\beta) + H(\beta)(\beta^* - \beta) \]
Estimasi parameter:
\[ \hat{\beta} \approx \beta^{(t)} - H^{-1}(\beta^{(t)}) U(\beta^{(t)}) \]
Diagnostik dilakukan untuk menilai apakah model GLM yang dibangun sudah sesuai dan memberikan hasil yang valid. Evaluasi dapat dilakukan dengan pendekatan formal maupun visual.
1. Statistik Devians
Devians mengukur perbandingan antara model saat ini dengan saturated model:
\[ D = 2 \sum_i \left[ y_i \log \left( \frac{y_i}{\hat{\mu}_i} \right) - (y_i - \hat{\mu}_i) \right] \]
2. Statistik Chi-Kuadrat Pearson
Digunakan untuk mengevaluasi apakah model lebih baik daripada tidak ada model sama sekali:
\[ X^2 = \sum \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i} \]
Jika statistik ini signifikan, model dianggap memberikan penjelasan yang lebih baik daripada model kosong.
Catatan Tambahan
3. Analisis Residual
Residual adalah selisih antara nilai observasi dengan prediksi
model.
- Dapat digunakan untuk mendeteksi penyimpangan sistematis
- Dapat diplot untuk mengevaluasi asumsi model
Sub-bab ini berfokus pada regresi logistik, yaitu salah satu bentuk GLM yang digunakan untuk model biner. Penjelasan mencakup fungsi model, estimasi parameter dengan Newton-Raphson, serta uji inferensial menggunakan Uji Wald.
1. Fungsi Model Logistik dan Log-Likelihood
Model logistik digunakan untuk memodelkan probabilitas dari suatu kejadian biner:
\[ \pi(x) = \frac{\exp(\beta_0 + \beta_1 x)}{1 + \exp(\beta_0 + \beta_1 x)} \]
Log-likelihood untuk \(n\) observasi:
\[ \ell(\beta) = \sum_{i=1}^n [y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i)] \]
2. Estimasi dengan Newton-Raphson
Model logistik diestimasi dengan memaksimalkan fungsi log-likelihood menggunakan metode Newton-Raphson. Model probabilitas:
\[ \pi_i = \frac{1}{1 + \exp(-x_i^T \beta)} \]
Turunan pertama (score function):
\[ U(\beta) = \frac{\partial \ell(\beta)}{\partial \beta} = \mathbf{X}^T (y - \pi) \]
Turunan kedua (Hessian matrix):
\[ H(\beta) = -\mathbf{X}^T \mathbf{W} \mathbf{X}, \quad \text{dengan } \mathbf{W} = \text{diag}(\pi_i (1 - \pi_i)) \]
Iterasi:
\[ \beta^{(t+1)} = \beta^{(t)} + (\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T (y - \pi^{(t)}) \]
Implementasi Newton-Raphson di R
# Simulasi Data
set.seed(1)
x <- rnorm(100)
beta_true <- c(-1, 2)
X <- cbind(1, x)
eta <- X %*% beta_true
p <- 1 / (1 + exp(-eta))
y <- rbinom(100, 1, p)
# Newton-Raphson
beta <- c(0, 0)
tol <- 1e-6
max_iter <- 100
for (i in 1:max_iter) {
eta <- X %*% beta
pi_hat <- 1 / (1 + exp(-eta))
W <- diag(as.numeric(pi_hat * (1 - pi_hat)))
z <- eta + solve(W) %*% (y - pi_hat)
beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
if (sum(abs(beta_new - beta)) < tol) break
beta <- beta_new
}
beta
## [,1]
## -1.516679
## x 2.155658
3. Uji Wald untuk Inferensi Parameter
Uji Wald digunakan untuk menguji signifikansi koefisien \(\beta_j\) dalam model regresi
logistik.
- Hipotesis:
\[ H_0: \beta_j = 0 \]
\[ H_1: \beta_j \ne 0 \]
- Asumsi: \[\hat{\beta}_j \sim N(\beta_j,
\text{Var}(\hat{\beta}_j))\]
- Statistik:
\[ Z = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \sim \mathcal{N}(0, 1) \] \[ W = Z^2 \sim \chi^2_1 \]
Simulasi dan Perhitungan Langkah-demi-Langkah
# Model regresi logistik dan uji Wald
set.seed(123)
n <- 100
x <- rnorm(n)
log_odds <- -0.5 + 1.2 * x
p <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, p)
data <- data.frame(x, y)
model <- glm(y ~ x, data = data, family = binomial)
summary(model)
##
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3097 0.2296 -1.349 0.177
## x 1.2663 0.3080 4.111 3.94e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 137.99 on 99 degrees of freedom
## Residual deviance: 114.76 on 98 degrees of freedom
## AIC: 118.76
##
## Number of Fisher Scoring iterations: 4
# Langkah perhitungan Wald
beta_hat <- coef(model)["x"]
se_beta <- summary(model)$coefficients["x", "Std. Error"]
Z <- beta_hat / se_beta
Wald_stat <- Z^2
p_value <- 1 - pchisq(Wald_stat, df = 1)
Interpretasi Hasil
Kesimpulan: uji Wald menguji apakah suatu prediktor berpengaruh secara statistik dalam model GLM, khususnya regresi logistik.
4. Uji Likelihood Ratio (Chi-Square)
Uji rasio likelihood, yaitu metode statistik untuk membandingkan kecocokan antara model penuh dan model null (tanpa prediktor) dalam regresi logistik. Uji ini dilakukan dengan membandingkan log-likelihood dari kedua model.
# Model null
model_null <- glm(y ~ 1, data = data, family = binomial)
# Likelihood ratio test
anova(model_null, model, test = "Chisq")
Jika nilai p yang dihasilkan < 0.05, maka model penuh secara signifikan lebih baik dalam menjelaskan variasi data dibandingkan model null.
Bab ini menjelaskan beberapa kriteria untuk mengevaluasi seberapa baik model regresi logistik yang dibangun.
1. Akaike Information Criterion (AIC)
AIC digunakan untuk mengevaluasi trade-off antara goodness-of-fit dan kompleksitas model. Semakin kecil nilai AIC, semakin baik model.
AIC(model)
## [1] 118.7598
2. Bayesian Information Criterion (BIC)
BIC mirip seperti AIC tetapi memberikan penalti lebih besar terhadap kompleksitas model. Cocok digunakan saat membandingkan banyak model dengan jumlah parameter berbeda.
BIC(model)
## [1] 123.9701
Catatan Penting:
Contoh Regresi Logistik Multiple
Bagian ini memberikan contoh penerapan regresi logistik dengan lebih
dari satu prediktor menggunakan dataset iris. Kita akan
memodelkan kemungkinan sebuah bunga merupakan spesies
setosa (1) atau bukan (0), berdasarkan panjang sepal dan
panjang petal.
1. Data dan Transformasi
# Menggunakan dataset iris
data(iris)
iris$target <- ifelse(iris$Species == "setosa", 1, 0)
summary(iris[, c("target", "Sepal.Length", "Petal.Length")])
## target Sepal.Length Petal.Length
## Min. :0.0000 Min. :4.300 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:5.100 1st Qu.:1.600
## Median :0.0000 Median :5.800 Median :4.350
## Mean :0.3333 Mean :5.843 Mean :3.758
## 3rd Qu.:1.0000 3rd Qu.:6.400 3rd Qu.:5.100
## Max. :1.0000 Max. :7.900 Max. :6.900
2. Model Regresi Logistik Fit
model_logit <- glm(target ~ Sepal.Length + Petal.Length, data = iris, family = binomial)
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(model_logit)
##
## Call:
## glm(formula = target ~ Sepal.Length + Petal.Length, family = binomial,
## data = iris)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -78.26 316177.72 0.000 1.000
## Sepal.Length 39.83 81738.89 0.000 1.000
## Petal.Length -48.60 43659.68 -0.001 0.999
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1.9095e+02 on 149 degrees of freedom
## Residual deviance: 5.0543e-09 on 147 degrees of freedom
## AIC: 6
##
## Number of Fisher Scoring iterations: 25
3. Uji Signifikansi Model (G² / Likelihood Ratio Test)
model_null <- glm(target ~ 1, data = iris, family = binomial)
anova(model_null, model_logit, test = "Chisq")
4. Statistik Wald untuk Koefisien Individu
library(broom)
## Warning: package 'broom' was built under R version 4.4.3
tidy(model_logit)
Output menunjukkan koefisien regresi, standar error, nilai statistik Wald, dan p-value untuk masing-masing prediktor.
5. Interpretasi Odds Ratio dan Confidence Interval
coefs <- summary(model_logit)$coefficients
OR <- exp(coefs[, "Estimate"])
CI_lower <- exp(coefs[, "Estimate"] - 1.96 * coefs[, "Std. Error"])
CI_upper <- exp(coefs[, "Estimate"] + 1.96 * coefs[, "Std. Error"])
odds_ratio_table <- cbind(OR, CI_lower, CI_upper)
round(odds_ratio_table, 3)
## OR CI_lower CI_upper
## (Intercept) 0.000000e+00 0 Inf
## Sepal.Length 1.993132e+17 0 Inf
## Petal.Length 0.000000e+00 0 Inf
6. Goodness of Fit (Hosmer-Lemeshow Test)
library(ResourceSelection)
## Warning: package 'ResourceSelection' was built under R version 4.4.3
## ResourceSelection 0.3-6 2023-06-27
hoslem.test(iris$target, fitted(model_logit), g = 3)
## Warning in hoslem.test(iris$target, fitted(model_logit), g = 3): The data did
## not allow for the requested number of bins.
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: iris$target, fitted(model_logit)
## X-squared = 2.5271e-09, df = 0, p-value < 2.2e-16
Output akan menunjukkan statistik uji dan p-value. Jika p-value besar, maka model dianggap memiliki kecocokan yang baik terhadap data.
7. Kesimpulan
Model regresi logistik menunjukkan bahwa variabel
Sepal.Length dan Petal.Length berkontribusi
dalam membedakan spesies setosa dari spesies lain. Jika
koefisien signifikan, maka kedua prediktor secara statistik membantu
dalam memprediksi status target. Odds ratio membantu dalam memahami
kekuatan dan arah pengaruh masing-masing prediktor terhadap peluang
menjadi setosa.
Sub-bab ini menjelaskan bagaimana menghitung statistik G² (likelihood ratio test) secara manual pada regresi logistik biner.
1. Rumus Umum
\[ G^2 = -2(\ell_0 - \ell_1) \]-
\(\ell_0\): log-likelihood dari model
null (tanpa prediktor)
- \(\ell_1\): log-likelihood dari model
penuh (dengan prediktor)
2. Fungsi Log-Likelihood
\[ \ell = \sum_{i=1}^n [y_i \log(p_i) + (1 - y_i) \log(1 - p_i)] \]
Contoh Data
# Data dummy
data <- data.frame(
y = c(1, 0, 1, 0),
x1 = c(2, 1, 3, 2)
)
Model Null (Intercept Saja)
# Probabilitas tetap
p_null <- mean(data$y) # 0.5
loglik_null <- sum(data$y * log(p_null) + (1 - data$y) * log(1 - p_null))
loglik_null
## [1] -2.772589
# [1] -2.772589
Model Full (Dengan Prediktor)
# Misal hasil estimasi model: beta0 = -1, beta1 = 1
eta <- -1 + 1 * data$x1
p_hat <- 1 / (1 + exp(-eta))
loglik_full <- sum(data$y * log(p_hat) + (1 - data$y) * log(1 - p_hat))
loglik_full
## [1] -2.446599
# [1] -2.4466
Hitung G²
G2 <- -2 * (loglik_null - loglik_full)
G2
## [1] 0.6519803
# [1] 0.6519803
Interpretasi
Bandingkan nilai \(G^2\) dengan distribusi \(\chi^2\) dengan derajat bebas (df) = jumlah prediktor:
pchisq(G2, df = 1, lower.tail = FALSE)
## [1] 0.4194056
# [1] 0.4194056
Jika \(p < 0.05\), maka model dengan prediktor secara signifikan lebih baik.
Persamaan logit: \[ \log\left(\frac{p_i}{1 - p_i}\right) = \beta_0 + \beta_1 x_{1i} + \dots + \beta_k x_{ki} \]
Fungsi likelihood: \[ L = \prod_{i=1}^n p_i^{y_i} (1 - p_i)^{1 - y_i} \]
Log-likelihood: \[ \ell = \sum_{i=1}^n [y_i \log(p_i) + (1 - y_i) \log(1 - p_i)] \]
Langkah-langkah Manual
Contoh Kasus (Ilustrasi Perhitungan Manual)
Data:
Model Null: \[ p = \frac{1 + 0 + 1 + 0}{4} = 0.5 \] \[ \ell_0 = \sum_{i=1}^4 [y_i \log(0.5) + (1 - y_i) \log(0.5)] = 4 \cdot \log(0.5) = -2.7726 \]
Model Full (misal \(\hat{\beta}_0 = -1, \hat{\beta}_1 = 1\)):
| x1 | y | \(\hat{p}\) | loglik component |
|---|---|---|---|
| 2 | 1 | 0.7311 | -0.3133 |
| 1 | 0 | 0.5000 | -0.6931 |
| 3 | 1 | 0.8808 | -0.1269 |
| 2 | 0 | 0.7311 | -1.3133 |
\[ \ell_1 = -0.3133 + (-0.6931) + (-0.1269) + (-1.3133) = -2.4466 \]
Hitung G²: \[ G^2 = -2(\ell_0 - \ell_1) = -2(-2.7726 + 2.4466) = 0.652 \]
Sub-bab ini menjelaskan penerapan model regresi Poisson, metode estimasi parameter menggunakan IRLS (Iteratively Reweighted Least Squares), serta uji statistik seperti Wald dan Likelihood Ratio Test.
Model ini digunakan untuk memodelkan data jumlah kejadian (count) yang mengikuti distribusi Poisson. Misalnya, jumlah pelanggan, kecelakaan, atau insiden per waktu/ruang.
Fungsi Distribusi Poisson:
\[ P(Y_i = y_i) = \frac{e^{-\lambda_i} \lambda_i^{y_i}}{y_i!} \]
Model:
\[ \log(\lambda_i) = \mathbf{x}_i^T \beta \]
Fungsi log-likelihood: \[ \ell(\beta) = \sum_{i=1}^{n} \left[y_i \log(\lambda_i) - \lambda_i - \log(y_i!) \right] \] Dengan: \[ \lambda_i = \exp(\mathbf{x}_i^T \beta) \]
Estimasi MLE dilakukan secara iteratif menggunakan pendekatan IRLS.
Langkah-Langkah IRLS:
Simulasi Estimasi Poisson di R (IRLS Manual)
set.seed(123)
n <- 100
x <- rnorm(n)
X <- cbind(1, x) # Tambahkan intercept
beta_true <- c(0.5, 0.8)
eta <- X %*% beta_true
lambda <- exp(eta)
y <- rpois(n, lambda)
data <- data.frame(x = x, y = y)
Iterasi IRLS Manual
beta <- c(0, 0)
tol <- 1e-6
max_iter <- 100
for (i in 1:max_iter) {
eta <- X %*% beta
lambda <- exp(eta)
W <- diag(as.numeric(lambda))
z <- eta + (y - lambda) / lambda
beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
if (sum(abs(beta_new - beta)) < tol) {
cat("Konvergen pada iterasi ke-", i, "\n")
break
}
beta <- beta_new
}
## Konvergen pada iterasi ke- 8
beta # hasil estimasi akhir
## [,1]
## 0.4494951
## x 0.8600048
Perbandingan dengan glm()
model_glm <- glm(y ~ x, family = poisson)
summary(model_glm)
##
## Call:
## glm(formula = y ~ x, family = poisson)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.44950 0.08872 5.066 4.05e-07 ***
## x 0.86000 0.07463 11.523 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 245.05 on 99 degrees of freedom
## Residual deviance: 106.78 on 98 degrees of freedom
## AIC: 325.76
##
## Number of Fisher Scoring iterations: 5
Uji Wald
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)
cat("Z:", wald_z,
"\nChi-Square:", wald_chisq,
"\np-value:", p_value)
## Z: 11.52289
## Chi-Square: 132.777
## p-value: 0
Uji Likelihood Ratio
model_null <- glm(y ~ 1, family = poisson, data = data)
anova(model_null, model_glm, test = "Chisq")
Evaluasi Model (AIC & BIC)
AIC(model_glm)
## [1] 325.7582
BIC(model_glm)
## [1] 330.9685
Kesimpulan:
glm().Regresi logistik merupakan salah satu metode analisis statistik yang digunakan untuk memodelkan hubungan antara satu variabel respons biner (dua kategori) dengan satu atau lebih variabel prediktor. Dalam penerapannya, prediktor yang digunakan dapat memiliki skala pengukuran berbeda, yaitu:
Nominal: Variabel yang tidak memiliki urutan logis antar kategorinya, hanya sebagai pembeda. Contoh: jenis kelamin (laki-laki, perempuan). Dalam analisis regresi logistik, data nominal diubah menjadi variabel dummy (misal: Female = 1, Male = 0).
Ordinal: Variabel yang memiliki urutan logis antar kategori, tetapi jarak antar kategori belum tentu sama. Contoh: tingkat pendidikan (SMA < Sarjana < Master < Doktor). Dalam analisis, variabel ordinal bisa diperlakukan:
Rasio: Variabel numerik kontinu yang memiliki nol absolut dan rasio bermakna. Contoh: pendapatan (dalam juta rupiah per bulan). Data ini dapat langsung digunakan dalam model tanpa transformasi khusus.
Dalam membangun model regresi logistik dengan prediktor campuran (nominal, ordinal, rasio), kita perlu menentukan perlakuan masing-masing tipe variabel:
Model regresi logistik secara umum:
\[ \log \left( \frac{p}{1-p} \right) = \beta_0 + \beta_1X_1 + \beta_2X_2 + \cdots + \beta_kX_k \]
Dimana:
- \(p\) adalah probabilitas kejadian.
- \(X_1, X_2, \ldots, X_k\) adalah variabel prediktor.
Setiap variabel dikodekan sesuai dengan tipenya.
library(tibble)
## Warning: package 'tibble' was built under R version 4.4.3
# Simulasi variabel
n <- 600
# Variabel baru: status perokok
smoker_status <- sample(c("NonSmoker", "Smoker"), n, replace = TRUE)
# Variabel baru: level aktivitas fisik
activity_level <- sample(c("Low", "Medium", "High"),
n, replace = TRUE, prob = c(0.3, 0.5, 0.2))
# Variabel baru: umur
age <- rnorm(n, mean = 35, sd = 10)
# Model logit dengan variabel baru
logit_p <- -1.8 + 0.6 * (smoker_status == "Smoker") +
0.9 * as.numeric(factor(activity_level, ordered = TRUE)) +
0.04 * age
# Konversi logit menjadi probabilitas
p <- 1 / (1 + exp(-logit_p))
# Set seed untuk reproduktifitas
set.seed(456)
# Simulasi data binomial
health_outcome <- rbinom(n, 1, p)
# Buat tabel data
df_simulasi <- tibble(health_outcome, smoker_status, activity_level, age)
# Lihat 6 data pertama
head(df_simulasi)
Interpretasi: Dataset berisi 600 observasi dengan variabel smoker_status, activity_level, age, dan health_outcome
# Grouping dan summarising
df_simulasi %>%
dplyr::group_by(health_outcome) %>%
dplyr::summarise(
Jumlah = dplyr::n(),
Rata2_Age = mean(age)
)
Interpretasi: Distribusi jumlah individu yang sehat dan tidak sehat berdasarkan simulasi, serta rata-rata umur pada masing-masing kelompok.
Dalam regresi logistik, perlakuan terhadap variabel ordinal menjadi penting untuk menentukan bagaimana variabel tersebut akan dimasukkan dalam model. Variabel ordinal memiliki sifat berurutan, namun jarak antar level tidak harus sama. Ada dua pendekatan utama untuk menangani variabel ordinal:
Dalam pendekatan ini, variabel ordinal diperlakukan sebagai variabel kategori biasa (nominal). Untuk setiap kategori, dibuat variabel dummy (0 atau 1). Jika ada \(k\) kategori, maka akan dibentuk \(k-1\) variabel dummy, dengan satu kategori sebagai referensi.
Misalnya, untuk variabel tingkat pendidikan dengan kategori: Rendah, Menengah, Tinggi, kita dapat membuat dua variabel dummy:
Kategori “Rendah” berfungsi sebagai baseline.
Keuntungan pendekatan ini: Fleksibel, tidak mengasumsikan hubungan linear antara kategori.
Kelemahannya: Mengabaikan informasi urutan kategori.
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.4.3
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:vcdExtra':
##
## summarise
## 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
# Mengubah variabel activity_level menjadi faktor berurutan
df_simulasi_nominal <- df_simulasi %>%
mutate(
activity_level = factor(activity_level, levels = c("Low", "Medium", "High"))
)
# Membuat model regresi logistik
model_nominal <- glm(health_outcome ~ smoker_status + activity_level + age,
data = df_simulasi_nominal, family = binomial)
# Menampilkan ringkasan model
summary(model_nominal)
##
## Call:
## glm(formula = health_outcome ~ smoker_status + activity_level +
## age, family = binomial, data = df_simulasi_nominal)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.57397 0.45404 1.264 0.20619
## smoker_statusSmoker 0.66283 0.23412 2.831 0.00464 **
## activity_levelMedium 0.83640 0.27937 2.994 0.00275 **
## activity_levelHigh -0.77575 0.27675 -2.803 0.00506 **
## age 0.01767 0.01165 1.516 0.12950
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 540.67 on 599 degrees of freedom
## Residual deviance: 495.72 on 595 degrees of freedom
## AIC: 505.72
##
## Number of Fisher Scoring iterations: 5
Keterangan: Model menggunakan activity_level sebagai variabel dummy, baseline “Low”. Setiap koefisien membandingkan kategori terhadap baseline.
Interpretasi Koefisien
Intercept (+0.12735)
Ini adalah log-odds dasar untuk individu NonSmoker dengan aktivitas fisik “Low” dan age = 0.
Tidak signifikan (p = 0.79381), artinya baseline ini tidak berbeda secara signifikan dari probabilitas 50%.
Odds Ratio = 1.14: Peluang health_outcome sedikit lebih besar daripada baseline, namun tidak signifikan.
smoker_statusSmoker (+0.78879)
Individu dengan status Smoker memiliki log-odds health_outcome lebih tinggi sebesar 0.789 dibandingkan NonSmoker.
p = 0.00126 (signifikan di level 1%), menunjukkan bahwa perokok memiliki peluang health_outcome lebih tinggi.
Odds Ratio = 2.20: Peluang health_outcome perokok sekitar 2.2 kali dibandingkan non-perokok.
activity_levelMedium (+0.87791)
Individu dengan aktivitas fisik Medium memiliki log-odds health_outcome lebih tinggi sebesar 0.878 dibandingkan Low.
p = 0.00307 (signifikan), menunjukkan aktivitas fisik Medium meningkatkan peluang health_outcome dibanding Low.
Odds Ratio = 2.40: Peluang health_outcome individu dengan aktivitas fisik Medium sekitar 2.4 kali dibanding aktivitas Low.
activity_levelHigh (-0.92189)
Individu dengan aktivitas fisik High memiliki log-odds health_outcome lebih rendah sebesar 0.922 dibandingkan Low.
p = 0.00143 (signifikan), menunjukkan bahwa aktivitas fisik High berhubungan dengan penurunan peluang health_outcome.
Odds Ratio = 0.40: Peluang health_outcome individu dengan aktivitas fisik High sekitar 40% dari peluang aktivitas Low.
age (+0.03146)
Setiap kenaikan 1 tahun usia meningkatkan log-odds health_outcome sebesar 0.0315.
p = 0.01143 (signifikan), artinya usia berhubungan positif dengan peluang health_outcome.
Odds Ratio = 1.032: Setiap tambahan 1 tahun usia meningkatkan peluang health_outcome sekitar 3.2%.
Interpretasi Goodness-of-Fit
Null deviance (524.27): Deviance model tanpa prediktor.
Residual deviance (467.39): Deviance model dengan prediktor.
AIC (477.39):
Signifikansi Model
Variabel smoker_status, activity_level (Medium dan High), dan age signifikan mempengaruhi peluang health_outcome.
Individu dengan status perokok dan aktivitas fisik Medium meningkatkan peluang health_outcome.
Aktivitas fisik High berhubungan dengan penurunan peluang health_outcome.
Kesimpulan Praktis
Status perokok dan aktivitas fisik berpengaruh signifikan terhadap peluang health_outcome.
Aktivitas fisik Medium meningkatkan peluang health_outcome, namun aktivitas fisik High justru menurunkannya.
Usia yang lebih tua sedikit meningkatkan peluang health_outcome.
Model cukup baik dalam memprediksi health_outcome berdasarkan variabel-variabel yang tersedia.
Alternatif lain adalah mengubah kategori ordinal menjadi angka berdasarkan rankingnya. Misalnya:
Kemudian angka-angka ini dimasukkan sebagai variabel numerik ke dalam model.
Keuntungan:
- Menjaga informasi urutan kategori.
- Lebih hemat parameter dibandingkan dummy.
Kelemahan:
- Mengasumsikan jarak antar level adalah sama.
Pilihan pendekatan tergantung pada konteks dan asumsi yang wajar berdasarkan data.
# Mengubah variabel activity_level menjadi skala numerik
df_simulasi_numeric <- df_simulasi %>%
mutate(
activity_level_numeric = case_when(
activity_level == "Low" ~ 1,
activity_level == "Medium" ~ 2,
activity_level == "High" ~ 3
)
)
# Membuat model regresi logistik dengan variabel numerik
model_numeric <- glm(health_outcome ~ smoker_status + activity_level_numeric + age,
data = df_simulasi_numeric, family = binomial)
# Menampilkan ringkasan model
summary(model_numeric)
##
## Call:
## glm(formula = health_outcome ~ smoker_status + activity_level_numeric +
## age, family = binomial, data = df_simulasi_numeric)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.45645 0.52465 2.776 0.00550 **
## smoker_statusSmoker 0.65794 0.22840 2.881 0.00397 **
## activity_level_numeric -0.38522 0.15852 -2.430 0.01510 *
## age 0.01784 0.01125 1.585 0.11292
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 540.67 on 599 degrees of freedom
## Residual deviance: 523.24 on 596 degrees of freedom
## AIC: 531.24
##
## Number of Fisher Scoring iterations: 4
Interpretasi: Model menggunakan
activity_level sebagai angka bertingkat 1–3. Koefisien
activity_level_numeric mengukur efek kenaikan satu tingkat
aktivitas fisik terhadap peluang health_outcome.
list(
AIC_Nominal = AIC(model_nominal),
AIC_Numeric = AIC(model_numeric)
)
## $AIC_Nominal
## [1] 505.7235
##
## $AIC_Numeric
## [1] 531.2398
Keterangan: Model dengan nilai AIC lebih kecil lebih
disukai. Membandingkan mana yang lebih baik: treat
activity_level sebagai nominal atau numeric.
Interpretasi Koefisien
Intercept (+1.53821)
Ini adalah log-odds dasar untuk individu NonSmoker dengan activity_level = 1 (Low) dan age = 0.
Signifikan (p = 0.00474), artinya baseline ini berbeda secara signifikan dari probabilitas 50%.
Odds Ratio = 4.66: Peluang health_outcome dasar sekitar 4.66 kali lebih besar dibanding baseline, dan signifikan.
smoker_statusSmoker (+0.76158)
Individu Smoker memiliki log-odds health_outcome lebih tinggi sebesar 0.762 dibandingkan NonSmoker.
p = 0.00134 (signifikan di level 1%), menunjukkan bahwa perokok memiliki peluang health_outcome yang lebih tinggi.
Odds Ratio = 2.14: Peluang health_outcome perokok sekitar 2.14 kali dibandingkan non-perokok.
activity_level_numeric (-0.53229)
Setiap kenaikan satu tingkat aktivitas fisik (Low -> Medium -> High) menurunkan log-odds health_outcome sebesar 0.532.
p = 0.00165 (signifikan), menunjukkan bahwa semakin tinggi tingkat aktivitas fisik, semakin menurun peluang health_outcome.
Odds Ratio = 0.59: Setiap kenaikan satu tingkat aktivitas fisik mengurangi peluang health_outcome sekitar 41%.
age (+0.02509)
Setiap kenaikan 1 tahun usia meningkatkan log-odds health_outcome sebesar 0.025.
p = 0.03446 (signifikan di level 5%), artinya usia berhubungan positif dengan peluang health_outcome.
Odds Ratio = 1.025: Setiap tambahan 1 tahun usia meningkatkan peluang health_outcome sekitar 2.5%.
Interpretasi Goodness-of-Fit
Null deviance (524.27): Deviance model tanpa prediktor.
Residual deviance (498.20): Deviance model dengan prediktor.
AIC (506.2):
Signifikansi Model
Variabel smoker_status, activity_level_numeric, dan age signifikan dalam mempengaruhi peluang health_outcome.
Smoker meningkatkan peluang health_outcome, sedangkan peningkatan aktivitas fisik menurunkan peluang health_outcome.
Usia yang lebih tinggi sedikit meningkatkan peluang health_outcome.
Kesimpulan Praktis
activity_level sebagai nominal atau numeric perlu
dipertimbangkan untuk optimasi.# Model null untuk baseline comparison
nullmod <- glm(health_outcome ~ 1, data = df_simulasi, family = binomial)
# Menghitung McFadden's R2
r2_nominal <- 1 - (logLik(model_nominal) / logLik(nullmod))
r2_numeric <- 1 - (logLik(model_numeric) / logLik(nullmod))
# Menampilkan hasil R2
list(
McFadden_R2_Nominal = r2_nominal,
McFadden_R2_Numeric = r2_numeric
)
## $McFadden_R2_Nominal
## 'log Lik.' 0.08313704 (df=5)
##
## $McFadden_R2_Numeric
## 'log Lik.' 0.03224425 (df=4)
Interpretasi: McFadden’s R² mengukur goodness-of-fit. Semakin besar nilainya, semakin baik model memprediksi dibandingkan model null.
library(ggplot2)
# Menambahkan prediksi ke data
sim_data_nominal <- df_simulasi_nominal %>%
mutate(predicted = predict(model_nominal, type = "response"))
sim_data_numeric <- df_simulasi_numeric %>%
mutate(predicted = predict(model_numeric, type = "response"))
# Plot untuk model nominal
sim_data_nominal %>%
ggplot(aes(x = age, y = predicted, color = smoker_status)) +
geom_point(alpha = 0.6) +
labs(
title = "Prediksi Probabilitas (Ordinal sebagai Nominal)",
x = "Usia",
y = "Prediksi Probabilitas Health Outcome"
) +
theme_minimal()
# Plot untuk model numeric
sim_data_numeric %>%
ggplot(aes(x = age, y = predicted, color = as.factor(activity_level_numeric))) +
geom_point(alpha = 0.6) +
labs(
title = "Prediksi Probabilitas (Ordinal sebagai Numeric)",
x = "Usia",
y = "Prediksi Probabilitas Health Outcome"
) +
theme_minimal()
Interpretasi: Visualisasi hubungan antara usia dan probabilitas health outcome berdasarkan status perokok dan tingkat aktivitas fisik, dengan dua pendekatan perlakuan ordinal.
library(knitr)
library(kableExtra)
library(broom)
# Ringkasan model nominal
tidy(model_nominal) %>%
kable(
format = ifelse(knitr::is_latex_output(), "latex", "html"),
booktabs = TRUE,
caption = "Ringkasan Koefisien Model dengan Motivasi sebagai Nominal"
) %>%
kable_styling(
latex_options = c("hold_position", "striped"),
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.5739655 | 0.4540415 | 1.264125 | 0.2061850 |
| smoker_statusSmoker | 0.6628313 | 0.2341202 | 2.831158 | 0.0046380 |
| activity_levelMedium | 0.8363976 | 0.2793725 | 2.993844 | 0.0027549 |
| activity_levelHigh | -0.7757468 | 0.2767532 | -2.803027 | 0.0050625 |
| age | 0.0176671 | 0.0116531 | 1.516090 | 0.1294967 |
# Ringkasan model numeric
tidy(model_numeric) %>%
kable(
format = ifelse(knitr::is_latex_output(), "latex", "html"),
booktabs = TRUE,
caption = "Ringkasan Koefisien Model dengan Motivasi sebagai Numeric"
) %>%
kable_styling(
latex_options = c("hold_position", "striped"),
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE
)
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 1.4564540 | 0.5246484 | 2.776057 | 0.0055023 |
| smoker_statusSmoker | 0.6579425 | 0.2283951 | 2.880720 | 0.0039677 |
| activity_level_numeric | -0.3852163 | 0.1585237 | -2.430023 | 0.0150979 |
| age | 0.0178394 | 0.0112538 | 1.585197 | 0.1129215 |
Interpretasi: Tabel memperlihatkan estimasi koefisien, standard error, nilai z, dan p-value. Koefisien positif meningkatkan log-odds health outcome, koefisien negatif menurunkan.
Kesimpulan - Smoker Status: - Individu dengan status perokok memiliki peluang health outcome lebih tinggi dibandingkan non-perokok.
Model dengan perlakuan ordinal sebagai numeric memiliki interpretasi lebih sederhana jika asumsi jarak antar tingkatannya konsisten.
Catatan: Semua hasil berbasis data simulasi, sehingga hasil bisa berbeda jika menggunakan data nyata.
Dalam analisis regresi logistik, memilih model yang tepat adalah tahap penting untuk memastikan bahwa model yang dihasilkan mampu menjelaskan data dengan baik dan memberikan hasil yang dapat diinterpretasikan. Pemilihan model tidak hanya mempertimbangkan kebaikan model dalam menyesuaikan data, tetapi juga memperhatikan kesederhanaan dan kemampuan generalisasi model tersebut.
Proses pemilihan model melibatkan serangkaian langkah yang bertujuan untuk menemukan model terbaik berdasarkan keseimbangan antara kompleksitas dan performa prediktif. Selain itu, evaluasi model juga diperlukan untuk menguji sejauh mana model dapat diterapkan pada data baru, sehingga hasil analisis tidak hanya berlaku pada data pelatihan.
Kriteria umum dalam pemilihan model mencakup:
- Goodness of fit: Seberapa baik model menjelaskan
variabilitas dalam data.
- Parsimony: Memilih model yang sederhana namun tetap
mampu menjelaskan hubungan antar variabel.
- Signifikansi statistik: Menguji kontribusi
masing-masing prediktor.
- Validasi model: Menilai kinerja model pada data yang
tidak digunakan saat proses fitting.
Dalam membangun model regresi logistik, terdapat dua pendekatan utama yang dapat digunakan:
Pendekatan confirmatory digunakan ketika peneliti telah memiliki
hipotesis atau model teoritis yang jelas sebelum analisis dilakukan.
Dalam pendekatan ini:
- Model dikembangkan berdasarkan teori atau penelitian sebelumnya.
- Variabel-variabel yang dimasukkan ke dalam model telah ditentukan
sebelumnya.
- Tujuannya adalah untuk menguji apakah data empiris mendukung model
yang dihipotesiskan.
Keunggulan pendekatan confirmatory adalah meningkatkan fokus analisis karena variabel yang diuji sudah dipilih secara teoritis. Namun, pendekatan ini memiliki kelemahan jika model awal tidak cukup baik untuk menggambarkan pola data yang sebenarnya.
Pendekatan exploratory lebih fleksibel dan digunakan ketika tidak ada model teoritis yang kuat atau ketika peneliti ingin mengeksplorasi data untuk menemukan pola yang belum diketahui. Dalam pendekatan ini, proses pemilihan variabel dapat dilakukan dengan beberapa metode, yaitu:
Forward Selection: Dimulai dari model kosong tanpa prediktor. Variabel dimasukkan satu per satu berdasarkan kontribusi terbesar terhadap model hingga tidak ada lagi variabel yang meningkatkan kualitas model secara signifikan.
Backward Elimination: Dimulai dari model penuh yang berisi semua variabel. Variabel dihapus satu per satu berdasarkan kontribusi terkecil, mengeluarkan variabel yang tidak signifikan hingga semua variabel yang tersisa relevan.
Stepwise Selection: Kombinasi dari forward selection dan backward elimination. Proses ini menambahkan dan/atau menghapus variabel pada setiap langkah untuk menemukan kombinasi prediktor terbaik.
Kelebihan pendekatan exploratory adalah memungkinkan penemuan hubungan baru antara variabel tanpa asumsi awal yang kuat. Namun, pendekatan ini rentan terhadap overfitting jika terlalu banyak variabel dipertahankan tanpa validasi silang atau kriteria seleksi yang ketat.
Tujuan Membangun Model
Tujuan utama dalam membangun model regresi logistik adalah menemukan model yang parsimonious, yaitu cukup sederhana namun tetap memiliki performa prediksi yang baik.
Pemilihan antara pendekatan Confirmatory dan
Exploratory sangat bergantung pada tujuan penelitian: -
Jika ingin menguji hipotesis tertentu berdasarkan teori yang ada,
pendekatan confirmatory menjadi pilihan.
- Jika tujuan utamanya adalah menemukan model terbaik dari data yang
tersedia, pendekatan exploratory lebih disarankan.
Dalam praktiknya, kedua pendekatan ini sering digunakan secara komplementer: teori digunakan sebagai dasar pengembangan model awal dan metode eksploratori digunakan untuk menyempurnakan model tersebut berdasarkan data empiris.
Dengan mempertimbangkan kedua pendekatan ini, peneliti dapat memilih strategi yang paling sesuai dengan tujuan penelitian dan karakteristik data yang dianalisis.
Simulasi Data
library(knitr)
library(dplyr)
library(ggplot2)
library(MASS)
## Warning: package 'MASS' was built under R version 4.4.3
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(caret)
## Warning: package 'caret' was built under R version 4.4.3
## Loading required package: lattice
##
## Attaching package: 'lattice'
## The following object is masked from 'package:gnm':
##
## barley
##
## Attaching package: 'caret'
## The following objects are masked from 'package:DescTools':
##
## MAE, RMSE
library(pROC)
## Warning: package 'pROC' was built under R version 4.4.3
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(DescTools)
set.seed(456)
n <- 250
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.4)
x3 <- rnorm(n)
lin_pred <- -0.8 + 1.5 * x1 - 0.5 * 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)
Pemilihan Model
model_full <- glm(y ~ x1 + x2 + x3, data = df, family = binomial)
summary(model_full)
##
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.7653 0.2292 -3.339 0.000841 ***
## x1 1.4113 0.2052 6.878 6.08e-12 ***
## x2 -0.5925 0.3226 -1.837 0.066241 .
## x3 0.6169 0.1751 3.522 0.000428 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 321.83 on 249 degrees of freedom
## Residual deviance: 238.95 on 246 degrees of freedom
## AIC: 246.95
##
## Number of Fisher Scoring iterations: 5
null_model <- glm(y ~ 1, data = df, family = binomial)
step_forward <- step(null_model, direction = "forward", scope = formula(model_full), trace = FALSE)
step_backward <- step(model_full, direction = "backward", trace = FALSE)
step_both <- step(null_model, direction = "both", scope = formula(model_full), trace = FALSE)
AIC(model_full, step_forward, step_backward, step_both)
pred_prob <- predict(step_both, type = "response")
roc_obj <- roc(df$y, pred_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj, main = "Kurva ROC", col = "blue")
auc(roc_obj)
## Area under the curve: 0.8267
PseudoR2(step_both, which = c("CoxSnell", "Nagelkerke", "McFadden"))
## CoxSnell Nagelkerke McFadden
## 0.2821649 0.3897376 0.2575264
pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
conf_matrix <- confusionMatrix(factor(pred_class), df$y, positive = "1")
conf_matrix
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 141 35
## 1 23 51
##
## Accuracy : 0.768
## 95% CI : (0.7107, 0.8189)
## No Information Rate : 0.656
## P-Value [Acc > NIR] : 8.103e-05
##
## Kappa : 0.4683
##
## Mcnemar's Test P-Value : 0.1486
##
## Sensitivity : 0.5930
## Specificity : 0.8598
## Pos Pred Value : 0.6892
## Neg Pred Value : 0.8011
## Prevalence : 0.3440
## Detection Rate : 0.2040
## Detection Prevalence : 0.2960
## Balanced Accuracy : 0.7264
##
## 'Positive' Class : 1
##
conf_matrix$byClass[c("Sensitivity", "Specificity")]
## Sensitivity Specificity
## 0.5930233 0.8597561
Model terbaik berdasarkan AIC adalah model model_full,
step_forward, step_backward, dan
step_both dengan nilai AIC yang sama yaitu 246.9477. ROC
menunjukkan nilai Area Under Curve (AUC) sebesar 0.83, yang
mengindikasikan performa model cukup baik dalam membedakan antara kelas
positif dan negatif. Pseudo-R² dan hasil klasifikasi juga menunjukkan
bahwa model memiliki akurasi dan kesesuaian yang cukup baik.
Secara keseluruhan, model yang dihasilkan cukup stabil dan akurat berdasarkan kombinasi metrik AIC, ROC AUC, dan klasifikasi.
Bagian ini ini bertujuan untuk memberikan panduan dalam membandingkan model regresi logistik dengan menggunakan berbagai ukuran evaluasi, yaitu: Deviance, AIC (Akaike Information Criterion), dan Likelihood-Ratio Test. Selain itu, dokumen ini juga membahas prinsip Parsimony yang menekankan pentingnya memilih model yang paling sederhana tanpa mengorbankan kualitas prediksi.
library(MASS)
library(broom)
library(DescTools)
set.seed(456)
Simulasi Data
# Simulasi Data
n <- 350
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.4)
x3 <- rnorm(n)
lin_pred <- -0.7 + 1.5 * x1 - 0.4 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)
head(data)
Pembuatan Model
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)
Perbandingan AIC dan Devianc
model_comp <- data.frame(
Model = c("Model 1", "Model 2", "Model 3"),
AIC = c(AIC(model1), AIC(model2), AIC(model3)),
Deviance = c(deviance(model1), deviance(model2), deviance(model3))
)
model_comp
anova(model1, model2, test = "LRT")
anova(model2, model3, test = "LRT")
Model dengan kompleksitas tinggi cenderung memiliki nilai AIC dan deviance yang lebih kecil. Namun, penting diingat bahwa model yang terlalu rumit bisa sulit diinterpretasikan. Oleh karena itu, model sederhana biasanya lebih disukai. Apabila penurunan AIC tidak memberikan peningkatan yang signifikan, sebaiknya pilih model yang lebih sederhana untuk menghindari overfitting.
Rumus dan Penjelasan
Rumus AIC
\[ \text{AIC} = -2 (\log L - k) = -2 \log L + 2k \]
Penjelasan: AIC atau Akaike Information Criterion digunakan untuk mengevaluasi model berdasarkan kombinasi kualitas kesesuaian model (melalui nilai log-likelihood) dan kompleksitas model (jumlah parameter \(k\)). Nilai AIC yang lebih rendah menunjukkan model yang lebih baik, karena AIC memberikan penalti terhadap model yang terlalu kompleks meskipun likelihood-nya tinggi.
Rumus Deviance
\[ D = -2 \left[ \log L(\text{model}) - \log L(\text{model saturasi}) \right] \]
Penjelasan: Deviance adalah ukuran yang menggambarkan seberapa jauh model saat ini dari model sempurna atau saturasi. Semakin kecil nilai deviance, semakin baik model dalam merepresentasikan data aktual, karena menunjukkan prediksi model semakin mendekati hasil pengamatan.
Rumus Likelihood-Ratio
\[ G^2 = -2 (\log L_0 - \log L_1) \]
Penjelasan: Statistik Likelihood Ratio digunakan untuk membandingkan dua model, yaitu model yang lebih sederhana dengan model yang lebih kompleks. Nilai \(G^2\) yang lebih kecil, atau p-value yang signifikan, menunjukkan bahwa model kompleks memberikan peningkatan kecocokan model yang signifikan dibandingkan model dasar.
# Prediksi probabilitas dan klasifikasi
pred_prob <- predict(model3, type = "response")
pred_class <- factor(ifelse(pred_prob >= 0.6, 1, 0))
# Membuat confusion matrix
conf_matrix <- confusionMatrix(pred_class, data$y, positive = "1")
conf_matrix
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 208 55
## 1 15 72
##
## Accuracy : 0.8
## 95% CI : (0.7542, 0.8406)
## No Information Rate : 0.6371
## P-Value [Acc > NIR] : 2.525e-11
##
## Kappa : 0.536
##
## Mcnemar's Test P-Value : 3.141e-06
##
## Sensitivity : 0.5669
## Specificity : 0.9327
## Pos Pred Value : 0.8276
## Neg Pred Value : 0.7909
## Prevalence : 0.3629
## Detection Rate : 0.2057
## Detection Prevalence : 0.2486
## Balanced Accuracy : 0.7498
##
## 'Positive' Class : 1
##
\[ \text{Sensitivity} = \frac{TP}{TP + FN} \]
\[ \text{Specificity} = \frac{TN}{TN + FP} \]
conf_matrix$byClass[c("Sensitivity", "Specificity")]
## Sensitivity Specificity
## 0.5669291 0.9327354
Kesimpulan
Kurva ROC adalah representasi visual yang digunakan untuk menilai kinerja model klasifikasi biner. Kurva ini menggambarkan hubungan antara True Positive Rate (Sensitivity) dan False Positive Rate (1 - Specificity) pada berbagai nilai threshold.
Dengan menggunakan package pROC kita dapat membuat kurva ROC:
library(pROC)
set.seed(456)
# Simulasi data
x1 <- rnorm(250)
x2 <- rbinom(250, 1, 0.4)
x3 <- rnorm(250)
lin_pred <- -0.8 + 1.3 * x1 - 0.6 * x2 + 0.5 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(250, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)
# Membuat model dan ROC
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
pred <- predict(model, type = "response")
roc_obj <- roc(data$y, pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj)
auc(roc_obj)
## Area under the curve: 0.7973
Dalam memilih threshold terbaik, kita bisa mengevaluasi sensitivitas dan spesifisitas pada berbagai cut-off.
# Menghitung sensitivitas dan spesifisitas untuk berbagai threshold
thresholds <- seq(0.15, 0.85, by = 0.1)
results <- data.frame(Threshold = thresholds)
results$Sensitivity <- sapply(thresholds, function(t) {
pred_class <- ifelse(pred >= t, 1, 0)
cm <- table(Pred = pred_class, Obs = data$y)
TP <- cm["1", "1"]
FN <- cm["0", "1"]
TP / (TP + FN)
})
results$Specificity <- sapply(thresholds, function(t) {
pred_class <- ifelse(pred >= t, 1, 0)
cm <- table(Pred = pred_class, Obs = data$y)
TN <- cm["0", "0"]
FP <- cm["1", "0"]
TN / (TN + FP)
})
print(results)
## Threshold Sensitivity Specificity
## 1 0.15 0.9375 0.4294118
## 2 0.25 0.7875 0.6235294
## 3 0.35 0.6875 0.7411765
## 4 0.45 0.5750 0.8529412
## 5 0.55 0.4000 0.9294118
## 6 0.65 0.2625 0.9470588
## 7 0.75 0.1250 0.9705882
## 8 0.85 0.0000 0.9941176
Cut-off terbaik dapat dipilih dengan kriteria:
Penjelasan Precision-Recall Curve
Kurva Precision-Recall (PR) digunakan sebagai metode evaluasi performa model klasifikasi, terutama efektif ketika menangani data dengan ketidakseimbangan kelas (class imbalance).
\[ \text{Precision} = \frac{TP}{TP + FP} \]
\[ \text{Recall} = \frac{TP}{TP + FN} \]
| Aspek | ROC Curve | Precision-Recall Curve |
|---|---|---|
| Fokus | Semua kelas | Hanya kelas positif |
| Kelebihan | Cocok untuk data seimbang | Efektif untuk data tidak seimbang |
| Sumbu Y | Sensitivitas (Recall) | Presisi |
| Sumbu X | 1 - Spesifisitas | Recall |
library(PRROC)
## Warning: package 'PRROC' was built under R version 4.4.3
## Loading required package: rlang
set.seed(456)
# Simulasi data
a1 <- rnorm(250)
a2 <- rbinom(250, 1, 0.4)
a3 <- rnorm(250)
lin_pred <- -0.5 + 1.2 * a1 - 0.5 * a2 + 0.7 * a3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(250, 1, p)
data <- data.frame(y = y, a1, a2, a3)
# Membuat model logistik
model <- glm(y ~ a1 + a2 + a3, data = data, family = binomial)
# Prediksi probabilitas
prob <- predict(model, type = "response")
# Membuat kurva Precision-Recall
pr <- pr.curve(scores.class0 = prob[data$y == 1],
scores.class1 = prob[data$y == 0],
curve = TRUE)
# Plot kurva
plot(pr)
Bagian ini menjelaskan dan menghitung pseudo R-squared dalam regresi logistik.
Simulasi Data
set.seed(456)
n <- 350
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.4)
x3 <- rnorm(n)
lin_pred <- -0.8 + 1.5 * x1 - 0.5 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)
Model Logistik dan Null Model
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
model_null <- glm(y ~ 1, data = data, family = binomial)
Likelihood dan Rumus
\[ R^2_{\text{Cox and Snell}} = 1 - \left( \frac{L_0}{L_M} \right)^{2/n} \]
\[ R^2_{\text{McFadden}} = 1 - \frac{\log L_M}{\log L_0} \]
Dengan:
Perhitungan Manual R-squared
logL0 <- logLik(model_null)
logLM <- logLik(model)
L0 <- exp(logL0)
LM <- exp(logLM)
n <- nobs(model)
cox_snell <- 1- (L0 / LM)^(2 / n)
mcfadden <- 1- (as.numeric(logLM) / as.numeric(logL0))
r2 <- data.frame(
R2_Cox_Snell = cox_snell,
R2_McFadden = mcfadden
)
r2
Perhitungan Otomatis dengan Package Tambahan
Menggunakan pscl
# Pilih mirror CRAN
options(repos = c(CRAN = "https://cran.r-project.org"))
# Cek dan install jika belum ada
if (!require(pscl)) install.packages("pscl")
## Loading required package: pscl
## Warning: package 'pscl' was built under R version 4.4.3
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
library(pscl)
# Hitung pseudo R-squared untuk model regresi logistik
pR2(model)
## fitting null model for pseudo-r2
## llh llhNull G2 McFadden r2ML r2CU
## -147.4803120 -224.3624173 153.7642105 0.3426693 0.3555296 0.4920566
Menggunakan rcompanion
if (!require(rcompanion)) install.packages("rcompanion"); library(rcompanion)
## Loading required package: rcompanion
## Warning: package 'rcompanion' was built under R version 4.4.3
nagelkerke(model)
## $Models
##
## Model: "glm, y ~ x1 + x2 + x3, binomial, data"
## Null: "glm, y ~ 1, binomial, data"
##
## $Pseudo.R.squared.for.model.vs.null
## Pseudo.R.squared
## McFadden 0.342669
## Cox and Snell (ML) 0.355530
## Nagelkerke (Cragg and Uhler) 0.492057
##
## $Likelihood.ratio.test
## Df.diff LogLik.diff Chisq p.value
## -3 -76.882 153.76 4.0615e-33
##
## $Number.of.observations
##
## Model: 350
## Null: 350
##
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
##
## $Warnings
## [1] "None"
Menggunakan DescTools
if (!require(DescTools)) install.packages("DescTools"); library(DescTools)
PseudoR2(model, which = "all")
## McFadden McFaddenAdj CoxSnell Nagelkerke AldrichNelson
## 0.3426693 0.3248410 0.3555296 0.4920566 0.3052305
## VeallZimmermann Efron McKelveyZavoina Tjur AIC
## 0.5433067 0.3927130 0.5676043 0.3921749 302.9606240
## BIC logLik logLik0 G2
## 318.3923566 -147.4803120 -224.3624173 153.7642105
Interpretasi
Gunakan kedua ukuran ini sebagai pelengkap untuk memberikan gambaran yang lebih komprehensif tentang performa model logistik.
Distribusi multinomial dapat dianggap sebagai bentuk generalisasi dari distribusi binomial untuk kasus dengan lebih dari dua kategori.
Misalkan \(X_1, X_2, \ldots, X_k\) menunjukkan jumlah kejadian dalam masing-masing dari \(k\) kategori, maka probabilitas distribusi multinomial didefinisikan sebagai:
\[ P(X_1 = x_1, \ldots, X_k = x_k) = \frac{n!}{x_1! x_2! \ldots x_k!} p_1^{x_1} p_2^{x_2} \ldots p_k^{x_k} \]
dengan ketentuan: \[ \sum_{i=1}^{k} x_i = n \quad \text{dan} \quad \sum_{i=1}^{k} p_i = 1.\]
Artinya, jumlah seluruh kejadian harus sama dengan \(n\) dan total probabilitas dari semua kategori harus 1.
Distribusi ini sangat berguna untuk memodelkan percobaan yang hasilnya dapat dikelompokkan ke dalam lebih dari dua kategori yang saling eksklusif.
Sebuah lembaga penelitian melakukan survei terhadap 15 orang untuk mengetahui pilihan moda transportasi favorit mereka: Mobil (M), Sepeda Motor (S), dan Kereta Api (K).
Hasil survei yang diperoleh:
- 6 orang memilih Mobil - 5 orang memilih Sepeda Motor - 4 orang memilih Kereta Api
Berdasarkan studi sebelumnya, probabilitas teoritis untuk tiap pilihan adalah:
- \(p_M = 0.4\)
- \(p_S = 0.35\)
- \(p_K = 0.25\)
Pertanyaannya: Berapakah peluang bahwa dari 15 orang, 6 orang memilih Mobil, 5 orang memilih Sepeda Motor, dan 4 orang memilih Kereta Api?
Untuk menentukan peluang ini, digunakan rumus distribusi multinomial:
\[ P(X_1 = x_1, \ldots, X_k = x_k) = \frac{n!}{x_1! x_2! \ldots x_k!} p_1^{x_1} p_2^{x_2} \ldots p_k^{x_k} \]
dengan:
- \(n = 15\), \(x_1 = 6\), \(x_2 = 5\), \(x_3 = 4\)
- \(p_1 = 0.4\), \(p_2 = 0.35\), \(p_3 = 0.25\)
Menggunakan rumus ini, kita bisa menghitung probabilitas dari kombinasi hasil yang diperoleh berdasarkan preferensi teoritis masing-masing moda transportasi.
Perhitungan Manual di R
n <- 15
x <- c(6, 5, 4)
p <- c(0.4, 0.35, 0.25)
# Hitung komponen-koefisien
faktorial_total <- factorial(n)
faktorial_x <- prod(factorial(x))
koefisien <- faktorial_total / faktorial_x
# Hitung peluang
peluang <- koefisien * prod(p^x)
peluang
## [1] 0.05299499
Hasil perhitungan menunjukkan bahwa peluang 6 orang memilih Mobil, 5 memilih Sepeda Motor, dan 4 memilih Kereta Api (dengan probabilitas masing-masing 0.4, 0.35, dan 0.25) adalah 0.05299499. Distribusi multinomial digunakan untuk menghitung peluang dalam percobaan yang melibatkan beberapa kategori hasil. Rumus dasarnya merupakan generalisasi dari distribusi binomial untuk lebih dari dua kategori.
Model ini berguna untuk menganalisis hubungan antara satu variabel respons kategori (dengan lebih dari dua kategori) dan satu atau lebih variabel independen.
Misalkan \(Y\) adalah variabel respons yang memiliki \(K\) kategori, dengan kategori \(K\) dipilih sebagai kategori referensi (baseline). Maka persamaan logit untuk kategori ke-\(j\) dapat dituliskan sebagai:
\[ \log \left( \frac{P(Y = j)}{P(Y = K)} \right) = \beta_{j0} + \beta_{j1}x_1 + \cdots + \beta_{jp}x_p \]
untuk \(j = 1, 2, \ldots, K-1\).
Artinya, model ini membandingkan probabilitas setiap kategori \(j\) terhadap kategori referensi \(K\), menggunakan kombinasi linier dari prediktor \(x_1, \ldots, x_p\).
Baseline-category logit model merupakan salah satu pendekatan dalam regresi logistik yang digunakan untuk menganalisis variabel respon kategori nominal dengan lebih dari dua kategori. Pada model ini, salah satu kategori dijadikan sebagai baseline (acuan), dan probabilitas kategori lain dibandingkan terhadap baseline tersebut menggunakan bentuk logit.
Rumus untuk fungsi logit kategori \(j\) adalah:
\[ \log \left( \frac{\pi_j}{\pi_c} \right), \quad j = 1, \ldots, c - 1 \]
dengan: - \(\pi_j\) adalah probabilitas respon pada kategori \(j\) - \(\pi_c\) adalah probabilitas respon pada kategori baseline \(c\)
Sehingga, terdapat \((c-1)\) fungsi logit yang dihasilkan.
Catatan: Pemilihan kategori baseline dapat dilakukan secara manual, namun secara default di R, kategori terakhir akan dijadikan baseline.
Model Regresi
Jika hanya ada satu prediktor \(x\), maka bentuk umum model logit yang digunakan adalah:
\[ \log \left( \frac{\pi_j}{\pi_c} \right) = \alpha_j + \beta_j x, \quad j = 1, \ldots, c-1 \]
Contoh Kasus: 3 Kategori Respon
Misalkan variabel respons \(Y\) memiliki tiga kategori: \(Y \in \{1, 2, 3\}\), dengan kategori ke-3 dipilih sebagai baseline. Maka model logit yang terbentuk adalah:
\[ \log \left( \frac{\pi_1}{\pi_3} \right) = \alpha_1 + \beta_1 x \]
\[ \log \left( \frac{\pi_2}{\pi_3} \right) = \alpha_2 + \beta_2 x \]
Dalam hal ini, terdapat dua persamaan logit: satu membandingkan kategori 1 terhadap kategori 3, dan satu lagi membandingkan kategori 2 terhadap kategori 3.
Relasi Antar Kategori
Untuk menghitung logit antara kategori 1 dan kategori 2 secara langsung, dapat dilakukan dengan:
\[ \log \left( \frac{\pi_1}{\pi_2} \right) = \log \left( \frac{\pi_1 / \pi_3}{\pi_2 / \pi_3} \right) = \log \left( \frac{\pi_1}{\pi_3} \right) - \log \left( \frac{\pi_2}{\pi_3} \right) \]
Sehingga:
\[ (\alpha_1 + \beta_1 x) - (\alpha_2 + \beta_2 x) = (\alpha_1 - \alpha_2) + (\beta_1 - \beta_2) x \]
Ringkasan Model Baseline-category Logit
Implementasi di R: Model ini dapat diterapkan menggunakan fungsi
multinom()dari packagennet. Pengaturan kategori baseline bisa dilakukan dengan fungsirelevel().
Proses estimasi parameter pada model ini dilakukan menggunakan pendekatan maximum likelihood, yang diselesaikan melalui algoritma iteratif seperti Newton-Raphson.
Fungsi log-likelihood yang digunakan adalah:
\[ \ell(\beta) = \sum_{i=1}^{n} \sum_{j=1}^{K} y_{ij} \log(\pi_{ij}) \]
dengan:
- \(\pi_{ij} = P(Y_i = j \mid x_i)\) adalah probabilitas individu ke-\(i\) berada pada kategori \(j\) berdasarkan prediktor \(x_i\). - \(y_{ij} = 1\) jika \(Y_i = j\), dan 0 jika tidak.
Estimasi ini mencari nilai parameter \(\beta\) yang memaksimalkan nilai fungsi log-likelihood tersebut.
Sebuah perusahaan retail ingin memahami faktor-faktor yang memengaruhi preferensi pelanggan dalam memilih jenis produk: Produk A, Produk B, atau Produk C.
Perusahaan ini melakukan survei terhadap 64 pelanggan dan mengumpulkan data sebagai berikut:
Tujuan dari penelitian ini adalah:
Mengetahui bagaimana pendapatan pelanggan berpengaruh terhadap preferensi dalam memilih jenis produk.
1. Input Data
customer_data <- data.frame(
income = c(25, 30, 28, 32, 31, 36, 38, 40, 42, 43,
45, 46, 47, 48, 50, 52, 54, 55, 57, 58,
60, 62, 64, 65, 67, 68, 70, 71, 73, 75,
77, 78, 80, 82, 83, 85, 87, 88, 90, 92,
95, 97, 98, 100, 102, 105, 107, 110, 112, 115,
118, 120, 122, 125, 128, 130, 132, 135, 138, 140,
143, 145, 148, 150),
purchase = factor(c("A", "B", "B", "A", "C", "A", "C", "B", "A", "C",
"A", "C", "C", "B", "A", "A", "B", "A", "B", "C",
"A", "C", "C", "B", "A", "C", "B", "A", "A", "C",
"B", "B", "A", "A", "C", "A", "A", "C", "B", "B",
"C", "C", "B", "C", "B", "C", "B", "C", "B", "C",
"A", "C", "B", "C", "B", "C", "C", "B", "B", "C",
"B", "A", "C", "B"),
levels = c("A", "B", "C")) # response order
)
2. Eksplorasi Data
summary(customer_data)
## income purchase
## Min. : 25.00 A:18
## 1st Qu.: 53.50 B:22
## Median : 79.00 C:24
## Mean : 82.64
## 3rd Qu.:110.50
## Max. :150.00
library(ggplot2)
ggplot(customer_data, aes(x = income, fill = purchase)) +
geom_histogram(binwidth = 5, position = "stack", color = "black") +
labs(title = "Distribusi Pendapatan berdasarkan Pilihan Produk",
x = "Pendapatan (juta)", y = "Jumlah Pelanggan") +
theme_minimal()
3. Multinomial Logistic Regression
library(nnet)
## Warning: package 'nnet' was built under R version 4.4.3
# Model multinomial (default: baseline = Produk A)
model_mlr <- multinom(purchase ~ income, data = customer_data)
## # weights: 9 (4 variable)
## initial value 70.311186
## final value 67.229494
## converged
summary(model_mlr)
## Call:
## multinom(formula = purchase ~ income, data = customer_data)
##
## Coefficients:
## (Intercept) income
## B -1.349683 0.01997742
## C -1.194741 0.01921248
##
## Std. Errors:
## (Intercept) income
## B 0.8365782 0.010133281
## C 0.8143547 0.009960154
##
## Residual Deviance: 134.459
## AIC: 142.459
4. Interpretasi Koefisien
z_values <- summary(model_mlr)$coefficients / summary(model_mlr)$standard.errors
p_values <- 2 * (1 - pnorm(abs(z_values)))
coef_table <- cbind(summary(model_mlr)$coefficients,
"z value" = round(z_values, 2),
"p value" = round(p_values, 4))
knitr::kable(coef_table, caption = "Koefisien Model Multinomial Logistic Regression")
| (Intercept) | income | (Intercept) | income | (Intercept) | income | |
|---|---|---|---|---|---|---|
| B | -1.349683 | 0.0199774 | -1.61 | 1.97 | 0.1067 | 0.0487 |
| C | -1.194741 | 0.0192125 | -1.47 | 1.93 | 0.1423 | 0.0537 |
Perhatikan ada 3 kategori pada variabel respons
purchase:
Maka model membangun dua model logit terhadap baseline (A), yaitu:
Logit untuk B terhadap A: \[ \log\left( \frac{P(B)}{P(A)} \right) = \beta_0^{(B)} + \beta_1^{(B)} \cdot \text{income} \]
Logit untuk C terhadap A: \[ \log\left( \frac{P(C)}{P(A)} \right) = \beta_0^{(C)} + \beta_1^{(C)} \cdot \text{income} \]
Jadi untuk setiap kategori selain referensi (A), kita dapatkan:
incomeBerikut adalah hasil estimasi model
multinom(purchase ~ income):
Koefisien Model Multinomial Logistic Regression
| Kategori | Intercept | Income | Z.Intercept | Z.Income | p.Intercept | p.Income |
|---|---|---|---|---|---|---|
| B | -1.349683 | 0.01997742 | -1.61 | 1.97 | 0.1067 | 0.0487 |
| C | -1.194741 | 0.01921248 | -1.47 | 1.93 | 0.1423 | 0.0537 |
Interpretasi
Untuk Kategori B: Koefisien income =
0.01997742, p-value = 0.0487 → signifikan.
Peningkatan pendapatan berasosiasi dengan peningkatan peluang memilih
Produk B dibanding Produk A.
Untuk Kategori C: Koefisien income =
0.01921248 → peluang memilih Produk C dibanding Produk A juga meningkat
seiring kenaikan pendapatan, namun p-value = 0.0537, sehingga marginal
signifikan.
Interpretasi dalam Odds Ratio
Kategori B (Produk B vs Produk A)
income = 0.01997742\[ \text{OR} = e^{0.01997742} \approx 1.0202 \]
Interpretasi:
Setiap peningkatan pendapatan 1 juta, peluang memilih Produk B dibanding
Produk A meningkat sekitar 2.02%, dan
signifikan secara statistik (p-value =
0.0487).
Kategori C (Produk C vs Produk A)
income = 0.01921248\[ \text{OR} = e^{0.01921248} \approx 1.0194 \]
Interpretasi:
Setiap peningkatan pendapatan 1 juta, peluang memilih Produk C dibanding
Produk A meningkat sekitar 1.94%, dengan p-value =
0.0537 (hampir signifikan).
Ringkasan
Ringkasan Interpretasi Odds Ratio
| Perbandingan | Koef.Income | Odds.Ratio | Interpretasi | Signifikansi |
|---|---|---|---|---|
| B vs A (Produk A) | 0.01997742 | 1.0202 | Peluang memilih B meningkat dengan income | Signifikan (p = 0.0487) |
| C vs A (Produk A) | 0.01921248 | 1.0194 | Peluang memilih C meningkat dengan income | Marginal (p = 0.0537) |
5. Prediksi dan Visualisasi
income_seq <- data.frame(income = seq(25, 150, by = 1))
predicted_probs <- predict(model_mlr, newdata = income_seq, type = "probs")
pred_df <- cbind(income_seq, predicted_probs)
pred_long <- tidyr::pivot_longer(pred_df, cols = -income, names_to = "purchase", values_to = "prob")
ggplot(pred_long, aes(x = income, y = prob, color = purchase)) +
geom_line(size = 1.2) +
labs(title = "Probabilitas Pilihan Produk berdasarkan Pendapatan",
x = "Pendapatan (juta)", y = "Probabilitas") +
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Penghitungan G² (Likelihood Ratio Test)
Definisi
G² (Deviance) dihitung sebagai: \[ G^2 = -2 \left( \log L_{\text{null}} - \log L_{\text{full}} \right) \]
McFadden’s Pseudo R² dihitung dengan rumus: \[ R^2 = 1 - \frac{\log L_{\text{full}}}{\log L_{\text{null}}} \]
Keterangan:
Kita bandingkan model null vs model penuh:
income# Model null (tanpa prediktor)
model_null <- multinom(purchase ~ 1, data = customer_data, trace = FALSE)
# Model penuh (dengan prediktor income)
model_full <- multinom(purchase ~ income, data = customer_data, trace = FALSE)
# Log-likelihood masing-masing
LL_null <- logLik(model_null)
LL_full <- logLik(model_full)
# G² = -2 (LL_null - LL_full)
G2 <- -2 * (as.numeric(LL_null) - as.numeric(LL_full))
# Derajat kebebasan = jumlah tambahan parameter
df <- attr(LL_full, "df") - attr(LL_null, "df")
# p-value
p_value <- 1 - pchisq(G2, df)
# Output
cat("Nilai G² =", round(G2, 4), "\n")
## Nilai G² = 5.2722
cat("Derajat kebebasan =", df, "\n")
## Derajat kebebasan = 2
cat("p-value =", round(p_value, 4), "\n")
## p-value = 0.0716
Sebuah perusahaan retail ingin mengetahui faktor-faktor yang memengaruhi preferensi pelanggan dalam memilih produk: Produk A, Produk B, atau Produk C.
Perusahaan ini melakukan survei terhadap sejumlah pelanggan dan mengumpulkan data sebagai berikut:
Tujuan dari penelitian ini adalah:
Mengetahui bagaimana pengaruh wilayah tempat tinggal dan kelompok usia terhadap preferensi pelanggan dalam memilih jenis produk.
1. Input Data
# Buat tabel sesuai dengan konteks baru
library(tibble)
library(dplyr)
library(nnet)
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.4.3
purchase_data <- tribble(
~Region, ~AgeGroup, ~Preference, ~Count,
"Urban", "Young", "A", 200,
"Urban", "Young", "B", 50,
"Urban", "Young", "C", 70,
"Urban", "Adult", "A", 180,
"Urban", "Adult", "B", 60,
"Urban", "Adult", "C", 90,
"Rural", "Young", "A", 90,
"Rural", "Young", "B", 30,
"Rural", "Young", "C", 40,
"Rural", "Adult", "A", 60,
"Rural", "Adult", "B", 20,
"Rural", "Adult", "C", 30
)
# Ubah ke format faktor
purchase_data <- purchase_data %>%
mutate(
Preference = factor(Preference, levels = c("A", "B", "C")),
AgeGroup = factor(AgeGroup),
Region = factor(Region)
)
2. Ekspansi Data (Untuk multinom())
# Perluas data sesuai Count
expanded_data <- purchase_data %>%
uncount(weights = Count)
head(expanded_data)
3. Model Multinomial Logistic Regression
# Model multinomial: Preference sebagai response
model_mlr <- multinom(Preference ~ Region + AgeGroup, data = expanded_data, trace = FALSE)
summary(model_mlr)
## Call:
## multinom(formula = Preference ~ Region + AgeGroup, data = expanded_data,
## trace = FALSE)
##
## Coefficients:
## (Intercept) RegionUrban AgeGroupYoung
## B -0.9776908 -0.1618435 -0.2031025
## C -0.5939113 -0.1323111 -0.2875300
##
## Std. Errors:
## (Intercept) RegionUrban AgeGroupYoung
## B 0.1946979 0.1969692 0.1814939
## C 0.1712427 0.1738306 0.1590161
##
## Residual Deviance: 1777.625
## AIC: 1789.625
4. Interpretasi Koefisien dan p-value
# Z dan p-value
z_values <- summary(model_mlr)$coefficients / summary(model_mlr)$standard.errors
p_values <- 2 * (1 - pnorm(abs(z_values)))
coef_table <- cbind(summary(model_mlr)$coefficients,
"z value" = round(z_values, 2),
"p value" = round(p_values, 4))
knitr::kable(coef_table, caption = "Koefisien Model Multinomial Logistic Regression") %>%
kable_styling(latex_options = "scale_down")
| (Intercept) | RegionUrban | AgeGroupYoung | (Intercept) | RegionUrban | AgeGroupYoung | (Intercept) | RegionUrban | AgeGroupYoung | |
|---|---|---|---|---|---|---|---|---|---|
| B | -0.9776908 | -0.1618435 | -0.2031025 | -5.02 | -0.82 | -1.12 | 0e+00 | 0.4113 | 0.2631 |
| C | -0.5939113 | -0.1323111 | -0.2875300 | -3.47 | -0.76 | -1.81 | 5e-04 | 0.4466 | 0.0706 |
5. Interpretasi Umum
Referensi kategori untuk Preference adalah “A” (karena urutan pertama).
Koefisien dalam model ini merepresentasikan perubahan log-odds dalam memilih “B” atau “C” dibandingkan “A”, berdasarkan variabel Region dan AgeGroup.
6. Contoh Prediksi
# Prediksi probabilitas
new_data <- expand.grid(
Region = levels(expanded_data$Region),
AgeGroup = levels(expanded_data$AgeGroup)
)
new_data$prob <- predict(model_mlr, newdata = new_data, type = "probs")
# Gabungkan hasil
cbind(new_data, new_data$prob)
Catatan Jika ingin menambahkan interaction (Age * AgeGroup), formula dapat diubah menjadi:
# Model dengan interaksi
model_mlr <- multinom(Preference ~ Region * AgeGroup, data = expanded_data)
## # weights: 15 (8 variable)
## initial value 1010.723306
## iter 10 value 888.532251
## final value 888.423667
## converged
Perhitungan G² dan McFadden Pseudo R²
# Model null: hanya intercept
model_null <- multinom(Preference ~ 1, data = expanded_data, trace = FALSE)
# Log-likelihood
ll_null <- logLik(model_null)
ll_full <- logLik(model_mlr)
# G2 (deviance)
G2 <- -2 * (as.numeric(ll_null) - as.numeric(ll_full))
# Derajat kebebasan
df <- attr(ll_full, "df") - attr(ll_null, "df")
# p-value
pval <- 1 - pchisq(G2, df)
# McFadden's pseudo R2
pseudo_r2 <- 1 - (as.numeric(ll_full) / as.numeric(ll_null))
# Tampilkan hasil
cat("Nilai G² (Deviance):", round(G2, 4), "\n")
## Nilai G² (Deviance): 5.1785
cat("Derajat Kebebasan:", df, "\n")
## Derajat Kebebasan: 6
cat("p-value:", round(pval, 4), "\n")
## p-value: 0.5211
cat("McFadden's Pseudo R²:", round(pseudo_r2, 4), "\n")
## McFadden's Pseudo R²: 0.0029
Regresi logistik ordinal digunakan saat variabel respon \(Y\) bersifat berurutan (ordinal), seperti misalnya tingkat kepuasan: Rendah, Sedang, dan Tinggi.
Berbeda dengan:
Model yang digunakan adalah Cumulative Logit Model dengan asumsi proportional odds:
\[ \log\left( \frac{P(Y \leq j)}{P(Y > j)} \right) = \alpha_j + \beta x \]
Dengan:
Jika terdapat \(c\) kategori, maka akan ada \(c-1\) model logit kumulatif.
Koefisien \(\beta\) menunjukkan bagaimana pengaruh \(x\) terhadap peluang untuk berada pada kategori yang lebih rendah atau sama.
Odds Ratio (OR) dihitung dengan:
\[ \text{OR} = e^{\beta} \]
Sebuah perusahaan ingin memahami faktor-faktor yang mempengaruhi kepuasan pelanggan terhadap layanan keanggotaan: Gold atau Silver.
Perusahaan ini melakukan survei terhadap pelanggan dan mengumpulkan data berikut:
Tujuan dari penelitian ini adalah:
Mengetahui bagaimana jenis keanggotaan dan jenis kelamin mempengaruhi tingkat kepuasan pelanggan terhadap layanan.
1. Input Data
# Data baru
library(tibble)
library(dplyr)
library(tidyr)
library(MASS)
customer_opinion <- tribble(
~Gender, ~Membership, ~Satisfaction, ~Count,
"Female", "Gold", "Very Satisfied", 35,
"Female", "Gold", "Satisfied", 95,
"Female", "Gold", "Neutral", 76,
"Female", "Gold", "Dissatisfied", 28,
"Female", "Gold", "Very Dissatisfied", 6,
"Female", "Silver", "Very Satisfied", 5,
"Female", "Silver", "Satisfied", 15,
"Female", "Silver", "Neutral", 22,
"Female", "Silver", "Dissatisfied", 73,
"Female", "Silver", "Very Dissatisfied", 30,
"Male", "Gold", "Very Satisfied", 30,
"Male", "Gold", "Satisfied", 70,
"Male", "Gold", "Neutral", 53,
"Male", "Gold", "Dissatisfied", 24,
"Male", "Gold", "Very Dissatisfied", 8,
"Male", "Silver", "Very Satisfied", 2,
"Male", "Silver", "Satisfied", 7,
"Male", "Silver", "Neutral", 19,
"Male", "Silver", "Dissatisfied", 67,
"Male", "Silver", "Very Dissatisfied", 35
)
# Pastikan urutan faktor ordinal sesuai spektrum tingkat kepuasan
customer_opinion <- customer_opinion %>%
mutate(
Satisfaction = factor(Satisfaction,
levels = c("Very Satisfied", "Satisfied", "Neutral",
"Dissatisfied", "Very Dissatisfied"),
ordered = TRUE),
Gender = factor(Gender),
Membership = factor(Membership)
)
2. Ekspansi Data
# Perluas berdasarkan Count
expanded_data <- uncount(customer_opinion, weights = Count)
head(expanded_data)
3. Fit Model: Ordinal Logistic Regression
# Model proportional odds menggunakan polr
model_ordinal <- polr(Satisfaction ~ Gender + Membership, data = expanded_data, Hess = TRUE)
summary(model_ordinal)
## Call:
## polr(formula = Satisfaction ~ Gender + Membership, data = expanded_data,
## Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## GenderMale 0.1684 0.139 1.211
## MembershipSilver 2.5391 0.172 14.760
##
## Intercepts:
## Value Std. Error t value
## Very Satisfied|Satisfied -1.5952 0.1417 -11.2540
## Satisfied|Neutral 0.2519 0.1118 2.2526
## Neutral|Dissatisfied 1.6436 0.1343 12.2410
## Dissatisfied|Very Dissatisfied 3.7425 0.1902 19.6775
##
## Residual Deviance: 1880.316
## AIC: 1892.316
4. Uji Signifikansi dan p-value
# Hitung z dan p-value
coefs <- coef(summary(model_ordinal))
p_values <- pnorm(abs(coefs[, "t value"]), lower.tail = FALSE) * 2
cbind(coefs, "p value" = round(p_values, 4)) %>%
knitr::kable(caption = "Koefisien dan p-value dari Model Regresi Logistik Ordinal")
| Value | Std. Error | t value | p value | |
|---|---|---|---|---|
| GenderMale | 0.1683786 | 0.1390172 | 1.211207 | 0.2258 |
| MembershipSilver | 2.5390932 | 0.1720235 | 14.760153 | 0.0000 |
| Very Satisfied|Satisfied | -1.5952275 | 0.1417475 | -11.254008 | 0.0000 |
| Satisfied|Neutral | 0.2519242 | 0.1118354 | 2.252635 | 0.0243 |
| Neutral|Dissatisfied | 1.6436104 | 0.1342704 | 12.241050 | 0.0000 |
| Dissatisfied|Very Dissatisfied | 3.7425214 | 0.1901933 | 19.677461 | 0.0000 |
5. Interpretasi Awal
Model ini mengasumsikan efek linier kumulatif dari prediktor terhadap tingkat kepuasan pelanggan yang bersifat ordinal.
6. Prediksi Probabilitas (Opsional)
new_data <- expand.grid(Gender = levels(expanded_data$Gender),
Membership = levels(expanded_data$Membership))
# Prediksi probabilitas
predict_probs <- predict(model_ordinal, newdata = new_data, type = "probs")
cbind(new_data, predict_probs)
Pada model cumulative logit, diasumsikan bahwa pengaruh prediktor \(x\) bersifat konstan di seluruh batas kategori (cutoff). Apabila asumsi ini tidak terpenuhi, disarankan untuk menggunakan pendekatan alternatif seperti generalized ordinal model.
Selain cumulative logit, terdapat beberapa model ordinal lain, antara lain:
Model-model ini biasanya dipilih ketika asumsi proportional odds pada cumulative logit tidak valid atau tidak sesuai dengan data yang dianalisis.
polr() dari paket MASS.Jika dibutuhkan validasi lebih lanjut, dapat dilakukan dengan uji devians atau likelihood ratio test untuk menguji ketepatan model.
Dalam regresi logistik ordinal, model yang sering dipakai adalah Cumulative Logit Model yang mengasumsikan Proportional Odds.
Asumsi ini juga dikenal dengan istilah parallel lines assumption, yang berarti:
\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j + \beta x \]
dengan \(j = 1, \ldots, c-1\), di mana: - Intercept (\(\alpha_j\)) berbeda untuk tiap batas kategori, - Koefisien \(\beta\) sama di seluruh kategori.
Visualisasi: Jika asumsi ini valid, maka kurva logit kumulatif terhadap prediktor akan memiliki kemiringan (slope) yang sama, berbeda hanya di titik potongnya (intercept).
Apa yang terjadi jika asumsi dilanggar? - Efek prediktor berbeda pada tiap batas kategori. - Model cumulative logit menjadi tidak sah. - Disarankan menggunakan model lain seperti: - Generalized Ordinal Logistic Regression - Partial Proportional Odds Model
Uji Asumsi Paralelisme Untuk mengecek validitas asumsi ini, dapat digunakan:
- Likelihood Ratio Test — membandingkan model proportional vs non-proportional. - Brant Test — tersedia dalam paket brant di R.
Ringkasan
Dalam statistika terapan, analisis data kategorik memiliki peran yang sangat penting karena banyak fenomena di dunia nyata berupa data dalam bentuk kategori, seperti jenis kelamin, tingkat pendidikan, preferensi konsumen, atau hasil diagnosis penyakit. Untuk mengolah data seperti ini, pendekatan yang banyak digunakan antara lain tabel kontingensi, model log-linear, dan regresi logistik. Setiap metode ini memiliki karakteristik serta kelebihan masing-masing, tergantung pada tujuan analisis dan pola data yang dianalisis.
Sebagai langkah awal eksplorasi, tabel kontingensi biasanya digunakan untuk menggambarkan hubungan antara dua atau lebih variabel kategorik. Misalnya, dalam studi tentang hubungan antara kebiasaan makan dan kejadian diabetes, tabel kontingensi dapat menampilkan frekuensi penderita diabetes berdasarkan jenis diet yang diikuti. Selain itu, tabel ini memudahkan perhitungan ukuran asosiasi seperti odds ratio, risk ratio, serta uji chi-square untuk menguji ada tidaknya hubungan antar variabel.
Apabila analisis memerlukan pengendalian terhadap banyak variabel sekaligus dan mempertimbangkan interaksinya, maka model log-linear menjadi alat yang sangat bermanfaat. Model ini merupakan bagian dari keluarga Generalized Linear Model (GLM) dan umumnya digunakan untuk menganalisis data berupa hitungan atau frekuensi dengan asumsi distribusi Poisson. Berbeda dari regresi logistik yang mengutamakan hubungan sebab akibat, model log-linear memperlakukan semua variabel secara setara, tanpa membedakan mana variabel dependen atau independen. Karena sifat simetris ini, model log-linear lebih sesuai untuk mengkaji struktur asosiasi dan independensi antar variabel.
Penyusunan model log-linear melibatkan pemodelan efek utama dari masing-masing variabel serta interaksi antar variabel tersebut. Misalnya, dalam studi yang meneliti hubungan antara tingkat pendidikan, status merokok, dan kejadian kanker paru-paru, model log-linear dapat digunakan untuk mengidentifikasi apakah interaksi tiga arah diperlukan untuk memahami pola hubungan di antara ketiga variabel tersebut. Untuk menentukan model yang paling tepat, metode likelihood ratio test kerap digunakan guna membandingkan model yang lebih sederhana dengan model yang lebih kompleks.
Sementara itu, regresi logistik banyak digunakan ketika satu variabel kategorik dijadikan sebagai variabel dependen, seperti kejadian stroke (ya/tidak), dan diprediksi oleh satu atau lebih variabel bebas, baik kategorik maupun numerik. Regresi logistik memodelkan log odds terjadinya suatu kejadian dan sangat berguna dalam penelitian observasional maupun eksperimental, terutama ketika outcome berbentuk biner. Model ini juga memiliki variasi lain, seperti regresi logistik multinomial untuk outcome dengan lebih dari dua kategori, atau regresi logistik ordinal ketika kategori memiliki urutan yang alami.
Secara keseluruhan, meskipun ketiga metode tersebut sama-sama menangani data kategorik, pendekatan dan fokus analisisnya berbeda. Tabel kontingensi berfungsi sebagai alat deskriptif awal, model log-linear digunakan untuk mengkaji struktur hubungan antar variabel secara simultan, sedangkan regresi logistik difokuskan pada pemodelan hubungan kausal dengan hasil akhir yang eksplisit.
Ringkasan
Dalam analisis data kategorik, pendekatan statistik yang umum meliputi:
Walaupun sama-sama diaplikasikan untuk data kategorik, masing-masing pendekatan memiliki interpretasi dan tujuan analisis yang berbeda.
Tabel Kontingensi
Tabel kontingensi digunakan untuk menyajikan frekuensi kombinasi antar kategori variabel.
Sebagai contoh, berikut adalah tabel 2x2:
plant_data <- matrix(c(45, 25, 55, 75), nrow = 2,
dimnames = list(Treatment = c("Fertilizer A", "Fertilizer B"),
Growth = c("Good", "Poor")))
plant_data
## Growth
## Treatment Good Poor
## Fertilizer A 45 55
## Fertilizer B 25 75
## Growth
## Treatment Good Poor
## Fertilizer A 45 25
## Fertilizer B 55 75
Tabel ini bersifat deskriptif dan tidak melibatkan pemodelan probabilitas.
Model Log-linear
Model log-linear memodelkan logaritma ekspektasi frekuensi sel dalam tabel kontingensi:
\[ \log(\mu_{ij}) = \mu + \lambda^{A}_i + \lambda^{B}_j + \lambda^{AB}_{ij} \]
Model log-linear dapat diterapkan dengan R sebagai berikut:
library(MASS)
loglm(~ Treatment * Growth, data = plant_data)
## Call:
## loglm(formula = ~Treatment * Growth, data = plant_data)
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 0 0 1
## Pearson 0 0 1
Model Regresi Logistik
Model regresi logistik biner digunakan untuk memodelkan probabilitas outcome kategorik:
\[ \log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1x \]
Contoh implementasi dalam R:
data_glm <- data.frame(
Growth = c(1, 0, 1, 0),
Treatment = factor(c("Fertilizer A", "Fertilizer A", "Fertilizer B", "Fertilizer B")),
Count = c(45, 25, 55, 75)
)
model_logit <- glm(Growth ~ Treatment, weights = Count, family = binomial, data = data_glm)
summary(model_logit)
##
## Call:
## glm(formula = Growth ~ Treatment, family = binomial, data = data_glm,
## weights = Count)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5878 0.2494 2.356 0.01845 *
## TreatmentFertilizer B -0.8979 0.3062 -2.933 0.00336 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 277.26 on 3 degrees of freedom
## Residual deviance: 268.38 on 2 degrees of freedom
## AIC: 272.38
##
## Number of Fisher Scoring iterations: 4
Perbandingan Tiga Pendekatan
| Aspek | Tabel Kontingensi | Model Log-linear | Regresi Logistik |
|---|---|---|---|
| Tujuan | Menyajikan frekuensi | Mengidentifikasi asosiasi | Memprediksi probabilitas |
| Variabel Dependen | Tidak ada | Tidak ada (hubungan simetris) | Ada (eksplisit) |
| Distribusi | Tidak diasumsikan | Poisson (frekuensi sel) | Binomial (probabilitas) |
| Bentuk Model | Tidak ada | GLM: log(\(\mu\)) ~ efek | GLM: logit(\(p\)) ~ prediktor |
| Cocok Untuk | Eksplorasi awal | Data tabel multi-dimensi (>2 variabel) | Analisis prediktif |
Tabel kontingensi digunakan untuk menampilkan frekuensi kombinasi kategori dari dua variabel atau lebih. Berikut contoh penggunaannya:
# Membuat tabel frekuensi untuk dua kategori
plant_data <- matrix(c(45, 25, 55, 75), nrow = 2,
dimnames = list(Treatment = c("Fertilizer A", "Fertilizer B"),
Growth = c("Good", "Poor")))
plant_data
## Growth
## Treatment Good Poor
## Fertilizer A 45 55
## Fertilizer B 25 75
Model log-linier untuk tabel \(I \times J\) dapat diformulasikan sebagai:
\[ \log(\mu_{ij}) = \lambda + \lambda^{T}_i + \lambda^{G}_j \]
Model saturated atau model penuh mencakup semua efek utama dan interaksi antara variabel:
Contoh penyusunan model untuk tabel 2x2:
# Import library
library(MASS)
# Membuat data
growth_data <- matrix(c(40, 60, 50, 50), nrow = 2, byrow = TRUE)
dimnames(growth_data) <- list(Treatment = c("Fertilizer A", "Fertilizer B"),
Growth = c("Good", "Poor"))
ftable(growth_data)
## Growth Good Poor
## Treatment
## Fertilizer A 40 60
## Fertilizer B 50 50
# Membuat model saturated
model_full <- loglm(~ Treatment * Growth, data = growth_data)
summary(model_full)
## Formula:
## ~Treatment * Growth
## attr(,"variables")
## list(Treatment, Growth)
## attr(,"factors")
## Treatment Growth Treatment:Growth
## Treatment 1 0 1
## Growth 0 1 1
## attr(,"term.labels")
## [1] "Treatment" "Growth" "Treatment:Growth"
## attr(,"order")
## [1] 1 1 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 0 0 1
## Pearson 0 0 1
Model saturated ini dapat dibangun menggunakan fungsi
loglm dari paket {MASS}.
Model tanpa interaksi antar variabel:
\[ \log(\mu_{ij}) = \mu + \lambda^{T}_i + \lambda^{G}_j \] Model ini menguji hipotesis bahwa variabel X dan Y saling independen.
model_indep <- loglm(~ Treatment + Growth, data = growth_data)
summary(model_indep)
## Formula:
## ~Treatment + Growth
## attr(,"variables")
## list(Treatment, Growth)
## attr(,"factors")
## Treatment Growth
## Treatment 1 0
## Growth 0 1
## attr(,"term.labels")
## [1] "Treatment" "Growth"
## 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 2.023756 1 0.1548557
## Pearson 2.020202 1 0.1552185
Odds ratio untuk tabel \(2 \times 2\) dihitung sebagai:
\[ OR = \frac{n_{11}n_{22}}{n_{12}n_{21}} \]
# Menghitung odds ratio
odds_ratio <- (growth_data[1,1] * growth_data[2,2]) / (growth_data[1,2] * growth_data[2,1])
odds_ratio
## [1] 0.6666667
Interpretasi:
Pada model saturated:
# Estimasi odds ratio dan log-odds
logOR <- log((growth_data[1,1] * growth_data[2,2]) / (growth_data[1,2] * growth_data[2,1]))
logOR
## [1] -0.4054651
Perbandingan model menggunakan deviance \(G^2\) atau likelihood ratio test:
anova(model_indep, model_full)
## LR tests for hierarchical log-linear models
##
## Model 1:
## ~Treatment + Growth
## Model 2:
## ~Treatment * Growth
##
## Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1 2.023756 1
## Model 2 0.000000 0 2.023756 1 0.15486
## Saturated 0.000000 0 0.000000 0 1.00000
Contoh studi hubungan antara tingkat pelayanan dan tingkat kepuasan pelanggan:
survey_data <- matrix(c(28, 172,
105, 595,
39, 341),
nrow = 3, byrow = TRUE,
dimnames = list(Service = c("Low", "Medium", "High"),
Satisfaction = c("Unsatisfied", "Satisfied")))
ftable(survey_data)
## Satisfaction Unsatisfied Satisfied
## Service
## Low 28 172
## Medium 105 595
## High 39 341
# Model log-linier untuk studi kasus
loglm(~ Service + Satisfaction, data = survey_data)
## Call:
## loglm(formula = ~Service + Satisfaction, data = survey_data)
##
## Statistics:
## X^2 df P(> X^2)
## Likelihood Ratio 5.019505 2 0.08128837
## Pearson 4.815504 2 0.09001741
Model log-linear adalah salah satu pendekatan untuk menganalisis hubungan antar dua atau lebih variabel kategori dalam sebuah tabel kontingensi. Dalam model ini, logaritma dari ekspektasi frekuensi sel (\(\mu_{ij}\)) diekspresikan sebagai penjumlahan efek variabel utama dan interaksinya. Untuk tabel 2x2, model ini ditulis sebagai:
\[ \log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^{AB}_{ij} \]
Sistem Persamaan Log-Linear
\[ \begin{aligned} \log(\mu_{11}) &= \lambda + \lambda^A_1 + \lambda^B_1 + \lambda^{AB}_{11} \\ \log(\mu_{12}) &= \lambda + \lambda^A_1 + \lambda^B_2 + \lambda^{AB}_{12} \\ \log(\mu_{21}) &= \lambda + \lambda^A_2 + \lambda^B_1 + \lambda^{AB}_{21} \\ \log(\mu_{22}) &= \lambda + \lambda^A_2 + \lambda^B_2 + \lambda^{AB}_{22} \end{aligned} \]
Kendala Sum-to-Zero
\[ \begin{aligned} \lambda^A_1 + \lambda^A_2 &= 0 \\ \lambda^B_1 + \lambda^B_2 &= 0 \\ \lambda^{AB}_{11} + \lambda^{AB}_{12} + \lambda^{AB}_{21} + \lambda^{AB}_{22} &= 0 \end{aligned} \]
Rumus Estimasi dengan Kendala Sum-to-Zero
\[ \begin{aligned} \lambda^A_1 = \frac{1}{2} \left( \log \mu_{11} + \log \mu_{12} - \log \mu_{21} - \log \mu_{22} \right) \\ \lambda^B_1 = \frac{1}{2} \left( \log \mu_{11} + \log \mu_{21} - \log \mu_{12} - \log \mu_{22} \right) \\ \lambda^{AB}_{12} = \frac{1}{4} \left( \log \mu_{12} - \log \mu_{11} - \log \mu_{22} + \log \mu_{21} \right) \end{aligned} \]
Diberikan data sebagai berikut:
| Mengkonsumsi Sayur | Sehat | Sakit |
|---|---|---|
| Rutin | 25 | 35 |
| Tidak Rutin | 15 | 30 |
Model log-linear untuk tabel \(2 \times 2\):
\[ \log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^{AB}_{ij} \]
Dengan asumsi kendala sum-to-zero:
\[ \sum_i \lambda^A_i = 0, \quad \sum_j \lambda^B_j = 0, \quad \sum_{i,j} \lambda^{AB}_{ij} = 0 \]
Misal:
Observasi:
Sistem persamaan:
\[ \begin{aligned} \log(\mu_{11}) = \lambda + \lambda^A_1 + \lambda^B_1 + \lambda^{AB}_{11} \\ \log(\mu_{12}) = \lambda + \lambda^A_1 + \lambda^B_2 + \lambda^{AB}_{12} \\ \log(\mu_{21}) = \lambda + \lambda^A_2 + \lambda^B_1 + \lambda^{AB}_{21} \\ \log(\mu_{22}) = \lambda + \lambda^A_2 + \lambda^B_2 + \lambda^{AB}_{22} \end{aligned} \]
Constraint sum-to-zero:
\[ \begin{aligned} \lambda^A_1 + \lambda^A_2 = 0 \\ \lambda^B_1 + \lambda^B_2 = 0 \\ \lambda^{AB}_{11} + \lambda^{AB}_{12} + \lambda^{AB}_{21} + \lambda^{AB}_{22} = 0 \end{aligned} \]
Langkah-langkah:
1. Hitung rata-rata log frekuensi sel:
\[ \lambda = \frac{1}{4} (\log 25 + \log 35 + \log 15 + \log 30) = 3.1918 \]
2. Efek utama A (Sayur):
\[ \lambda^A_1 = \frac{1}{2} \left( (\log 25 + \log 35) - (\log 15 + \log 30) \right) = \frac{1}{2} (6.5511 - 5.7038) = 0.4237 \]
\[ \lambda^A_2 = -0.4237 \]
3. Efek utama B (Status):
\[ \lambda^B_1 = \frac{1}{2} \left( (\log 25 + \log 15) - (\log 35 + \log 30) \right) = \frac{1}{2} (5.5215 - 6.4513) = -0.4649 \]
\[ \lambda^B_2 = 0.4649 \]
4. Efek interaksi:
\[ \lambda^{AB}_{11} = \frac{1}{4} (\log 25 - \log 35 - \log 15 + \log 30) = \frac{1}{4} (3.2189 - 3.5553 - 2.7081 + 3.4012) = 0.0899 \]
\[ \lambda^{AB}_{12} = -0.0899, \quad \lambda^{AB}_{21} = -0.0899, \quad \lambda^{AB}_{22} = 0.0899 \]
Ringkasan Parameter
Odds Ratio (\(OR\)) dihitung sebagai:
\[ OR = \frac{n_{11} n_{22}}{n_{12} n_{21}} = \frac{25 \times 30}{35 \times 15} = \frac{750}{525} = 1.4286 \]
Log odds ratio:
\[ \log(OR) = \log(1.4286) = 0.3550 \]
Standard error (SE):
\[ SE = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}} = \sqrt{\frac{1}{25} + \frac{1}{35} + \frac{1}{15} + \frac{1}{30}} = \sqrt{0.1686} = 0.4106 \]
95% Confidence Interval untuk log(OR):
\[ \log(OR) \pm 1.96 \times SE = 0.3550 \pm 1.96 \times 0.4106 \]
\[ = (0.3550 - 0.805), (0.3550 + 0.805) = (-0.4500, 1.1600) \]
Back-transform untuk mendapatkan CI dari OR:
\[ \text{Lower} = \exp(-0.4500) = 0.6376 \]
\[ \text{Upper} = \exp(1.1600) = 3.1903 \]
Jadi, \(OR \approx 1.43\) (95% CI: 0.64 – 3.19)
# Data 2x2
tabel <- matrix(c(25, 35, 15, 30), nrow = 2, byrow = TRUE)
colnames(tabel) <- c("Sehat", "Sakit")
rownames(tabel) <- c("Rutin", "Tidak Rutin")
tabel
## Sehat Sakit
## Rutin 25 35
## Tidak Rutin 15 30
# Ubah tabel menjadi data frame
survey_data <- as.data.frame(as.table(tabel))
colnames(survey_data) <- c("Sayur", "Status", "Freq")
survey_data
# Model tanpa interaksi
fit_no_inter <- glm(Freq ~ Sayur + Status, family = poisson, data = survey_data)
summary(fit_no_inter)
##
## Call:
## glm(formula = Freq ~ Sayur + Status, family = poisson, data = survey_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1293 0.1793 17.454 <2e-16 ***
## SayurTidak Rutin -0.2877 0.1972 -1.459 0.1446
## StatusSakit 0.4855 0.2010 2.416 0.0157 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8.92165 on 3 degrees of freedom
## Residual deviance: 0.76151 on 1 degrees of freedom
## AIC: 27.025
##
## Number of Fisher Scoring iterations: 4
# Model dengan interaksi
fit_inter <- glm(Freq ~ Sayur * Status, family = poisson, data = survey_data)
summary(fit_inter)
##
## Call:
## glm(formula = Freq ~ Sayur * Status, family = poisson, data = survey_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.2189 0.2000 16.094 <2e-16 ***
## SayurTidak Rutin -0.5108 0.3266 -1.564 0.118
## StatusSakit 0.3365 0.2619 1.285 0.199
## SayurTidak Rutin:StatusSakit 0.3567 0.4106 0.869 0.385
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8.9216e+00 on 3 degrees of freedom
## Residual deviance: 1.0658e-14 on 0 degrees of freedom
## AIC: 28.263
##
## Number of Fisher Scoring iterations: 3
Nilai \(\log(1.4286) = 0.3550\) dari odds ratio sejalan dengan estimasi efek interaksi dari output model di R.
Sebuah studi dilakukan untuk mengevaluasi hubungan antara Aktivitas Fisik (Rendah/Tinggi) dan Kategori Berat Badan (Kurus, Ideal, Gemuk).
| Aktivitas Fisik | Kurus | Ideal | Gemuk |
|---|---|---|---|
| Rendah | 15 | 22 | 12 |
| Tinggi | 18 | 25 | 8 |
Bentuk umum model log-linear untuk tabel 2x3 dengan kendala sum-to-zero:
\[ \log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^{AB}_{ij} \]
dengan:
Secara eksplisit:
\[ \log(\mu_{ij}) = \lambda + \lambda^A_1 (\text{Rendah}) + \lambda^A_2 (\text{Tinggi}) + \lambda^B_1 (\text{Kurus}) + \lambda^B_2 (\text{Ideal}) + \lambda^B_3 (\text{Gemuk}) + \lambda^{AB}_{ij} (\text{interaksi}) \]
# Membuat data frame dari tabel
tabel2x3 <- matrix(c(15, 22, 12, 18, 25, 8), nrow = 2, byrow = TRUE)
colnames(tabel2x3) <- c("Kurus", "Ideal", "Gemuk")
rownames(tabel2x3) <- c("Rendah", "Tinggi")
tabel2x3
## Kurus Ideal Gemuk
## Rendah 15 22 12
## Tinggi 18 25 8
# Konversi menjadi data frame
data2x3 <- as.data.frame(as.table(tabel2x3))
colnames(data2x3) <- c("Aktivitas", "BB", "Freq")
data2x3
# Model tanpa interaksi
fit_no_inter <- glm(Freq ~ Aktivitas + BB, family = poisson, data = data2x3)
summary(fit_no_inter)
##
## Call:
## glm(formula = Freq ~ Aktivitas + BB, family = poisson, data = data2x3)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.78316 0.20177 13.794 <2e-16 ***
## AktivitasTinggi 0.04001 0.20004 0.200 0.8415
## BBIdeal 0.35364 0.22711 1.557 0.1194
## BBGemuk -0.50078 0.28338 -1.767 0.0772 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 12.4712 on 5 degrees of freedom
## Residual deviance: 1.2301 on 2 degrees of freedom
## AIC: 36.799
##
## Number of Fisher Scoring iterations: 4
# Model dengan interaksi
fit_inter <- glm(Freq ~ Aktivitas * BB, family = poisson, data = data2x3)
summary(fit_inter)
##
## Call:
## glm(formula = Freq ~ Aktivitas * BB, family = poisson, data = data2x3)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.70805 0.25820 10.488 <2e-16 ***
## AktivitasTinggi 0.18232 0.34960 0.522 0.602
## BBIdeal 0.38299 0.33485 1.144 0.253
## BBGemuk -0.22314 0.38730 -0.576 0.565
## AktivitasTinggi:BBIdeal -0.05449 0.45572 -0.120 0.905
## AktivitasTinggi:BBGemuk -0.58779 0.57494 -1.022 0.307
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1.2471e+01 on 5 degrees of freedom
## Residual deviance: 1.5543e-15 on 0 degrees of freedom
## AIC: 39.569
##
## Number of Fisher Scoring iterations: 3
Model tanpa interaksi
Model dengan interaksi
Contoh Interpretasi:
Pada bagian ini kita akan melanjutkan pembahasan mengenai model log-linear, dengan fokus pada model yang melibatkan tiga variabel kategorik. Model ini memungkinkan adanya interaksi antara tiga variabel, yang membuat analisis menjadi lebih kompleks.
Model log-linear untuk tiga variabel kategorik (X, Y, Z) dapat dibangun dengan berbagai tingkat interaksi tergantung pada kompleksitas model yang ingin dibuat. Berikut adalah beberapa variasi model yang umum digunakan:
1. Model Saturated
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk} \]
Model saturated memasukkan semua kemungkinan efek interaksi, termasuk interaksi tiga arah (X, Y, Z).
2. Model Homogen
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]
Model ini hanya mempertimbangkan interaksi dua arah antar variabel dan tidak melibatkan interaksi tiga arah.
3. Model Conditional
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} \]
Model ini memuat interaksi antara X dengan Y dan X dengan Z.
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk} \]
Memuat interaksi antara Y dengan X dan Y dengan Z.
4. Model Joint Independence
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} \]
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{YZ}_{jk} \]
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} \]
5. Model Tanpa Interaksi
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k \]
Model ini hanya memasukkan efek utama tanpa mempertimbangkan interaksi antar variabel.
Dalam model log-linear tiga arah, pengujian interaksi dilakukan untuk menentukan ada atau tidaknya hubungan antar variabel.
Prosedurnya dilakukan secara hierarkis, mulai dari interaksi tiga arah hingga dua arah.
1. Pengujian Interaksi Tiga Arah (XYZ)
2. Pengujian Interaksi Dua Arah (XY, XZ, YZ)
Setiap tahap pengujian digunakan untuk menilai kecocokan model terhadap data yang diamati dan untuk menentukan struktur interaksi yang paling sesuai.
Menggunakan R untuk menyelesaikan studi kasus berikut:
Tabel berikut merupakan data hasil survei tahun 2000 mengenai jenis pekerjaan, tingkat pendidikan, dan sikap terhadap penggunaan teknologi baru. Susun dan interpretasikan model log-linear paling sederhana yang cocok untuk data ini. Jelaskan langkah-langkah yang diambil untuk memilih model terbaik serta asosiasi mana saja yang teridentifikasi. Tunjukkan pula prediksi nilai dari model tersebut.
| Pekerjaan | Pendidikan | Setuju | Tidak Setuju | Total |
|---|---|---|---|---|
| Pegawai | Rendah | 140 | 30 | 170 |
| Pegawai | Tinggi | 110 | 50 | 160 |
| Pegawai Total | 250 | 80 | 330 | |
| Profesional | Rendah | 160 | 40 | 200 |
| Profesional | Tinggi | 150 | 60 | 210 |
| Profesional Total | 310 | 100 | 410 | |
| Wirausaha | Rendah | 130 | 60 | 190 |
| Wirausaha | Tinggi | 120 | 70 | 190 |
| Wirausaha Total | 250 | 130 | 380 |
Keterangan:
Package yang Digunakan
library("epitools")
library("DescTools")
library("lawstat")
## Warning: package 'lawstat' was built under R version 4.4.3
Input Data
# Input data sesuai tabel praktikum
z.job <- factor(rep(c("Pegawai", "Profesional", "Wirausaha"), each = 4))
x.educ <- factor(rep(c("Rendah", "Tinggi"), each = 2, times = 3))
y.tech <- factor(rep(c("Setuju", "TidakSetuju"), times = 6))
counts <- c(140, 30, 110, 50, 160, 40, 150, 60, 130, 60, 120, 70)
# Membuat data frame
data <- data.frame(
Pekerjaan = z.job,
Pendidikan = x.educ,
Sikap = y.tech,
Frekuensi = counts
)
data
Membentuk Tabel Kontingensi 3 Arah
table3d <- xtabs(Frekuensi ~ Pekerjaan + Pendidikan + Sikap, data = data)
ftable(table3d)
## Sikap Setuju TidakSetuju
## Pekerjaan Pendidikan
## Pegawai Rendah 140 30
## Tinggi 110 50
## Profesional Rendah 160 40
## Tinggi 150 60
## Wirausaha Rendah 130 60
## Tinggi 120 70
Analisis Log-Linear: Tahap Pemodelan
Selanjutnya kita akan memodelkan tabel ini menggunakan berbagai model log-linear dan membandingkan kecocokan model untuk menemukan model paling sederhana yang cukup menggambarkan data.
# Mengatur referensi kategori
x.educ <- relevel(x.educ, ref = "Tinggi")
y.tech <- relevel(y.tech, ref = "TidakSetuju")
z.job <- relevel(z.job, ref = "Wirausaha")
Model Saturated
Model log-linear saturated mencakup semua interaksi hingga tiga arah:
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk} \]
# Fitting model saturated
model_saturated <- glm(counts ~ x.educ * y.tech * z.job,
family = poisson(link = "log"))
summary(model_saturated)
##
## Call:
## glm(formula = counts ~ x.educ * y.tech * z.job, family = poisson(link = "log"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2485 0.1195 35.545 < 2e-16
## x.educRendah -0.1542 0.1759 -0.876 0.380927
## y.techSetuju 0.5390 0.1504 3.584 0.000339
## z.jobPegawai -0.3365 0.1852 -1.817 0.069193
## z.jobProfesional -0.1542 0.1759 -0.876 0.380927
## x.educRendah:y.techSetuju 0.2342 0.2167 1.081 0.279917
## x.educRendah:z.jobPegawai -0.3567 0.2903 -1.229 0.219238
## x.educRendah:z.jobProfesional -0.2513 0.2695 -0.933 0.351030
## y.techSetuju:z.jobPegawai 0.2495 0.2274 1.097 0.272632
## y.techSetuju:z.jobProfesional 0.3773 0.2144 1.760 0.078399
## x.educRendah:y.techSetuju:z.jobPegawai 0.5178 0.3414 1.517 0.129333
## x.educRendah:y.techSetuju:z.jobProfesional 0.2358 0.3187 0.740 0.459334
##
## (Intercept) ***
## x.educRendah
## y.techSetuju ***
## z.jobPegawai .
## z.jobProfesional
## x.educRendah:y.techSetuju
## x.educRendah:z.jobPegawai
## x.educRendah:z.jobProfesional
## y.techSetuju:z.jobPegawai
## y.techSetuju:z.jobProfesional .
## x.educRendah:y.techSetuju:z.jobPegawai
## x.educRendah:y.techSetuju:z.jobProfesional
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2.6635e+02 on 11 degrees of freedom
## Residual deviance: 3.8636e-14 on 0 degrees of freedom
## AIC: 98.905
##
## Number of Fisher Scoring iterations: 3
# Melihat eksponensial dari koefisien
exp(model_saturated$coefficients)
## (Intercept)
## 70.0000000
## x.educRendah
## 0.8571429
## y.techSetuju
## 1.7142857
## z.jobPegawai
## 0.7142857
## z.jobProfesional
## 0.8571429
## x.educRendah:y.techSetuju
## 1.2638889
## x.educRendah:z.jobPegawai
## 0.7000000
## x.educRendah:z.jobProfesional
## 0.7777778
## y.techSetuju:z.jobPegawai
## 1.2833333
## y.techSetuju:z.jobProfesional
## 1.4583333
## x.educRendah:y.techSetuju:z.jobPegawai
## 1.6783217
## x.educRendah:y.techSetuju:z.jobProfesional
## 1.2659341
Model saturated ini membentuk hubungan penuh antara tingkat pendidikan, sikap terhadap teknologi, dan jenis pekerjaan. Model ini mempertimbangkan semua efek utama, interaksi dua arah, dan interaksi tiga arah.
library(knitr)
library(kableExtra)
# Contoh tabel
tabel <- data.frame(
Parameter = c("(Intercept)", "x.educRendah", "y.techSetuju", "z.jobPegawai", "z.jobProfesional",
"x.educRendah:y.techSetuju", "x.educRendah:z.jobPegawai", "x.educRendah:z.jobProfesional",
"y.techSetuju:z.jobPegawai", "y.techSetuju:z.jobProfesional",
"x.educRendah:y.techSetuju:z.jobPegawai", "x.educRendah:y.techSetuju:z.jobProfesional"),
Estimate = c(4.25, -0.15, 0.54, -0.34, -0.15, 0.23, -0.36, -0.25, 0.25, 0.38, 0.52, 0.24),
StdError = c(0.12, 0.18, 0.15, 0.19, 0.18, 0.22, 0.29, 0.27, 0.23, 0.21, 0.34, 0.32),
zvalue = c(35.55, -0.88, 3.58, -1.82, -0.88, 1.08, -1.23, -0.93, 1.10, 1.76, 1.52, 0.74),
pvalue = c("<2e-16 ***", "0.381", "0.00034 ***", "0.069 .", "0.381", "0.280", "0.219", "0.351", "0.273", "0.078 .", "0.129", "0.459")
)
kable(tabel, format = ifelse(knitr::is_latex_output(), "latex", "html"), booktabs = TRUE, linesep = "") %>%
kable_styling(latex_options = c("striped", "scale_down"))
| Parameter | Estimate | StdError | zvalue | pvalue |
|---|---|---|---|---|
| (Intercept) | 4.25 | 0.12 | 35.55 | <2e-16 *** |
| x.educRendah | -0.15 | 0.18 | -0.88 | 0.381 |
| y.techSetuju | 0.54 | 0.15 | 3.58 | 0.00034 *** |
| z.jobPegawai | -0.34 | 0.19 | -1.82 | 0.069 . |
| z.jobProfesional | -0.15 | 0.18 | -0.88 | 0.381 |
| x.educRendah:y.techSetuju | 0.23 | 0.22 | 1.08 | 0.280 |
| x.educRendah:z.jobPegawai | -0.36 | 0.29 | -1.23 | 0.219 |
| x.educRendah:z.jobProfesional | -0.25 | 0.27 | -0.93 | 0.351 |
| y.techSetuju:z.jobPegawai | 0.25 | 0.23 | 1.10 | 0.273 |
| y.techSetuju:z.jobProfesional | 0.38 | 0.21 | 1.76 | 0.078 . |
| x.educRendah:y.techSetuju:z.jobPegawai | 0.52 | 0.34 | 1.52 | 0.129 |
| x.educRendah:y.techSetuju:z.jobProfesional | 0.24 | 0.32 | 0.74 | 0.459 |
Catatan Interpretasi
Model log-linear homogenous hanya mempertimbangkan efek utama dan semua interaksi dua arah, tanpa melibatkan interaksi tiga arah. Model ini secara matematis dinyatakan sebagai:
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]
# Model Homogenous
model_homogenous <- glm(counts ~ x.educ + y.tech + z.job +
x.educ*y.tech + x.educ*z.job + y.tech*z.job,
family = poisson(link = "log"))
summary(model_homogenous)
##
## Call:
## glm(formula = counts ~ x.educ + y.tech + z.job + x.educ * y.tech +
## x.educ * z.job + y.tech * z.job, family = poisson(link = "log"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.31290 0.10528 40.966 < 2e-16 ***
## x.educRendah -0.29934 0.13679 -2.188 0.028649 *
## y.techSetuju 0.43502 0.12567 3.462 0.000537 ***
## z.jobPegawai -0.49224 0.15627 -3.150 0.001633 **
## z.jobProfesional -0.22341 0.14572 -1.533 0.125237
## x.educRendah:y.techSetuju 0.45415 0.13573 3.346 0.000820 ***
## x.educRendah:z.jobPegawai 0.01575 0.15190 0.104 0.917440
## x.educRendah:z.jobProfesional -0.09403 0.14385 -0.654 0.513319
## y.techSetuju:z.jobPegawai 0.48373 0.16878 2.866 0.004156 **
## y.techSetuju:z.jobProfesional 0.48807 0.15874 3.075 0.002108 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 266.3481 on 11 degrees of freedom
## Residual deviance: 2.3235 on 2 degrees of freedom
## AIC: 97.229
##
## Number of Fisher Scoring iterations: 4
Pengujian ini menggunakan nilai residual deviance dari model saturated dan model homogenous.
# Selisih deviance antar model
deviance_model <- model_homogenous$deviance - model_saturated$deviance
deviance_model
## [1] 2.323527
# Derajat bebas = db model homogenous - db model saturated
derajat_bebas <- model_homogenous$df.residual - model_saturated$df.residual
derajat_bebas
## [1] 2
chi_tabel <- qchisq(1 - 0.05, df = derajat_bebas)
chi_tabel
## [1] 5.991465
keputusan <- ifelse(deviance_model <= chi_tabel, "Terima H0", "Tolak H0")
keputusan
## [1] "Terima H0"
Interpretasi:
Pada tingkat signifikansi 5%, tidak ditemukan bukti yang cukup untuk menolak \(H_0\). Dengan demikian, tidak terdapat bukti adanya interaksi tiga arah antara pendidikan, sikap terhadap teknologi, dan pekerjaan.
Catatan:
Rangkuman
Pengujian untuk Menentukan Adanya Interaksi Tiga Arah (Model Saturated vs Model Homogenous)
\[ \Delta\text{Deviance} = \text{Deviance model homogenous} - \text{Deviance model saturated} \]
\[ = 2.324 - 0.00 = 2.324 \]
\[ \text{df} = \text{db model homogenous} - \text{db model saturated} = 2 \]
\[ \Delta\text{Deviance} > \chi^2_{0.05,2} = 5.991 \]
Catatan Tambahan: Perhitungan Derajat Bebas dan Selisih Deviance
Ketika menghitung selisih deviance, model yang lebih sederhana adalah model homogenous. Derajat bebas dihitung berdasarkan:
Periksa output R untuk memastikan jumlah parameter yang digunakan, termasuk intercept.
Model log-linear conditional pada X mempertimbangkan efek utama dan interaksi dua arah antara X dengan Y serta X dengan Z, tanpa memperhitungkan interaksi antara Y dan Z maupun interaksi tiga arah.
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} \]
# Conditional Association on X
model_conditional_X <- glm(counts ~ x.educ + y.tech + z.job +
x.educ*y.tech + x.educ*z.job,
family = poisson(link = "log"))
summary(model_conditional_X)
##
## Call:
## glm(formula = counts ~ x.educ + y.tech + z.job + x.educ * y.tech +
## x.educ * z.job, family = poisson(link = "log"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.11204 0.09504 43.266 < 2e-16 ***
## x.educRendah -0.32542 0.14214 -2.290 0.022050 *
## y.techSetuju 0.74721 0.09048 8.258 < 2e-16 ***
## z.jobPegawai -0.17185 0.10730 -1.602 0.109245
## z.jobProfesional 0.10008 0.10013 1.000 0.317512
## x.educRendah:y.techSetuju 0.44904 0.13493 3.328 0.000875 ***
## x.educRendah:z.jobPegawai 0.06062 0.15053 0.403 0.687135
## x.educRendah:z.jobProfesional -0.04879 0.14244 -0.343 0.731945
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 266.348 on 11 degrees of freedom
## Residual deviance: 14.477 on 4 degrees of freedom
## AIC: 105.38
##
## Number of Fisher Scoring iterations: 4
Pengujian Selisih Deviance (Conditional on X vs Homogenous) Langkah-langkah uji hipotesis menggunakan residual deviance:
# Deviasi antara model conditional X dan homogen
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance
Deviance.model
## [1] 12.15379
Menghitung Derajat Bebas
derajat.bebas <- (4 - 2)
derajat.bebas
## [1] 2
Menghitung Nilai Chi-Square Tabel
chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 5.991465
Keputusan Uji
keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
keputusan
## [1] "Tolak"
Interpretasi
Karena nilai \(\Delta\)Deviance = 2.13 lebih kecil dibandingkan nilai kritis chi-square tabel sebesar 5.991 (dengan df = 2, \(\alpha\) = 0.05), maka keputusan uji adalah Terima.
Kesimpulan
Pada tingkat signifikansi 5%, tidak terdapat cukup bukti untuk menolak \(H_0\). Dengan demikian, tidak terdapat interaksi signifikan antara sikap terhadap teknologi (Y) dan pekerjaan (Z). Model yang terbentuk hanya mencakup interaksi dua arah dengan variabel X (conditional on X) tanpa interaksi antara Y dan Z, yang berarti parameter \(\lambda^{YZ}_{jk}\) tidak signifikan.
Pengujian Ada Tidaknya Interaksi Antara Y dan Z (Homogenous Model vs Conditional on X)
Hipotesis
\(H_0\): \(\lambda^{YZ}_{jk} = 0\) (Tidak ada interaksi antara sikap terhadap teknologi (Y) dan pekerjaan (Z))
\(H_1\): \(\lambda^{YZ}_{jk} \neq 0\) (Ada interaksi antara sikap terhadap teknologi (Y) dan pekerjaan (Z))
Tingkat Signifikansi
Statistik Uji
\[ \Delta \text{Deviance} = \text{Deviance model conditional on X} - \text{Deviance model homogenous} \]
\[ = 14.477 - 2.324 = 12.153 \]
\[ \text{db} = \text{db model conditional on X} - \text{db model homogenous} = 4 - 2 = 2 \]
\[ \Delta \text{Deviance} > \chi^2_{0.05,2} = 5.991 \]
Model log-linear conditional pada Y mempertimbangkan efek utama dan interaksi dua arah antara X dengan Y serta Y dengan Z, tanpa memperhitungkan interaksi antara X dan Z maupun interaksi tiga arah.
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk} \]
# Conditional Association on Y
model_conditional_Y <- glm(counts ~ x.educ + y.tech + z.job +
x.educ*y.tech + y.tech*z.job,
family = poisson(link = "log"))
summary(model_conditional_Y)
##
## Call:
## glm(formula = counts ~ x.educ + y.tech + z.job + x.educ * y.tech +
## y.tech * z.job, family = poisson(link = "log"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.3239 0.1001 43.192 < 2e-16 ***
## x.educRendah -0.3254 0.1151 -2.827 0.004694 **
## y.techSetuju 0.4407 0.1242 3.549 0.000387 ***
## z.jobPegawai -0.4855 0.1421 -3.417 0.000634 ***
## z.jobProfesional -0.2624 0.1330 -1.972 0.048555 *
## x.educRendah:y.techSetuju 0.4490 0.1349 3.328 0.000875 ***
## y.techSetuju:z.jobPegawai 0.4855 0.1679 2.892 0.003833 **
## y.techSetuju:z.jobProfesional 0.4775 0.1579 3.025 0.002488 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 266.3481 on 11 degrees of freedom
## Residual deviance: 2.9937 on 4 degrees of freedom
## AIC: 93.899
##
## Number of Fisher Scoring iterations: 3
# Deviasi antara model conditional Y dan model homogenous
deviance_model <- model_conditional_Y$deviance - model_homogenous$deviance
deviance_model
## [1] 0.6701972
Menghitung Derajat Bebas
# Menghitung derajat bebas
derajat_bebas <- model_conditional_Y$df.residual - model_homogenous$df.residual
derajat_bebas
## [1] 2
Menghitung Nilai Chi-Square Tabel
# Menghitung nilai chi-square tabel dengan alpha = 0.05
chi_tabel <- qchisq(1 - 0.05, df = derajat_bebas)
chi_tabel
## [1] 5.991465
Keputusan Uji
# Keputusan Uji Hipotesis
keputusan <- ifelse(deviance_model <= chi_tabel, "Terima", "Tolak")
keputusan
## [1] "Terima"
Interpretasi
Karena nilai \(\Delta\)Deviance sebesar 0.670 lebih kecil dibandingkan nilai kritis chi-square tabel 5.991 (df = 2, \(\alpha\) = 0.05), maka keputusan uji adalah Terima.
Kesimpulan
Pada taraf signifikansi 5%, tidak terdapat cukup bukti untuk menolak \(H_0\). Artinya, tidak terdapat interaksi signifikan antara tingkat pendidikan (X) dan pekerjaan (Z) dalam model ini. Dengan demikian, model tanpa parameter interaksi \(\lambda^{XZ}_{ik}\) sudah memadai untuk menjelaskan data.
Pengujian Ada Tidaknya Interaksi Antara X dan Z (Homogenous Model vs Conditional Association on Y)
Hipotesis
\(H_0\): \(\lambda^{XZ}_{ik} = 0\) (Tidak ada interaksi antara tingkat pendidikan (X) dan pekerjaan (Z))
\(H_1\): \(\lambda^{XZ}_{ik} \neq 0\) (Ada interaksi antara tingkat pendidikan (X) dan pekerjaan (Z))
Tingkat Signifikansi
Statistik Uji
\[ \Delta \text{Deviance} = \text{Deviance model conditional on Y} - \text{Deviance model homogenous} \]
\[ = 2.994 - 2.324 = 0.670 \]
\[ \text{db} = \text{db model conditional on Y} - \text{db model homogenous} = 4 - 2 = 2 \]
Daerah Penolakan
\[ \Delta \text{Deviance} > \chi^2_{0.05,2} = 5.991 \]
Keputusan
Kesimpulan
Model log-linear conditional pada variabel Z mempertimbangkan efek utama dan interaksi dua arah antara X dengan Z dan Y dengan Z, tanpa memasukkan interaksi X dengan Y ataupun interaksi tiga arah.
Secara matematis, model ini dapat dituliskan:
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]
# Conditional Association on Z
model_conditional_Z <- glm(counts ~ x.educ + y.tech + z.job +
x.educ*z.job + y.tech*z.job,
family = poisson(link = "log"))
summary(model_conditional_Z)
##
## Call:
## glm(formula = counts ~ x.educ + y.tech + z.job + x.educ * z.job +
## y.tech * z.job, family = poisson(link = "log"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.174e+00 1.016e-01 41.084 < 2e-16 ***
## x.educRendah 5.433e-16 1.026e-01 0.000 1.00000
## y.techSetuju 6.539e-01 1.081e-01 6.048 1.47e-09 ***
## z.jobPegawai -5.163e-01 1.614e-01 -3.199 0.00138 **
## z.jobProfesional -2.383e-01 1.505e-01 -1.583 0.11336
## x.educRendah:z.jobPegawai 6.062e-02 1.505e-01 0.403 0.68713
## x.educRendah:z.jobProfesional -4.879e-02 1.424e-01 -0.343 0.73195
## y.techSetuju:z.jobPegawai 4.855e-01 1.679e-01 2.892 0.00383 **
## y.techSetuju:z.jobProfesional 4.775e-01 1.579e-01 3.025 0.00249 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 266.348 on 11 degrees of freedom
## Residual deviance: 13.635 on 3 degrees of freedom
## AIC: 106.54
##
## Number of Fisher Scoring iterations: 4
# Deviasi antara model conditional Z dan model homogen
deviance_model <- model_conditional_Z$deviance - model_homogenous$deviance
deviance_model
## [1] 11.31143
Menghitung Derajat Bebas
# Derajat bebas
derajat_bebas <- model_conditional_Z$df.residual - model_homogenous$df.residual
derajat_bebas
## [1] 1
Menghitung Nilai Chi-Square Tabel
# Nilai Chi-Square tabel
tabel_chi <- qchisq(1 - 0.05, df = derajat_bebas)
tabel_chi
## [1] 3.841459
Keputusan Uji
# Keputusan berdasarkan perbandingan deviance dan chi-square tabel
keputusan <- ifelse(deviance_model <= tabel_chi, "Terima", "Tolak")
keputusan
## [1] "Tolak"
Interpretasi
Karena nilai \(\Delta\)Deviance sebesar 27.9309 jauh lebih besar dari nilai kritis chi-square tabel 3.841 (df = 1, \(\alpha\) = 0.05), maka keputusan uji adalah Tolak.
Kesimpulan
Pada taraf nyata 5%, terdapat cukup bukti untuk menolak \(H_0\). Ini berarti, ada interaksi yang signifikan antara tingkat pendidikan (X) dan sikap terhadap teknologi (Y).
Pengujian Ada Tidaknya Interaksi Antara X dan Y (Homogen vs Conditional Association on Z)
Hipotesis
Tingkat Signifikansi
Statistik Uji
\[ \Delta \text{Deviance} = \text{Deviance model conditional on Z} - \text{Deviance model homogen} \]
\[ = 29.729 - 1.798 = 27.931 \]
\[ \text{db} = \text{db model conditional on Z} - \text{db model homogen} = 3 - 2 = 1 \]
Daerah Penolakan
\[ \Delta \text{Deviance} > \chi^2_{0.05,1} = 3.841 \]
Keputusan
Kesimpulan
| Model | Parameter | Deviance | Jumlah Parameter | df | AIC |
|---|---|---|---|---|---|
| Saturated | \(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk}\) | 0.00 | 12 | 0 | 100.14 |
| Homogen | \(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk}\) | 1.798 | 10 | 2 | 97.93 |
| Conditional on X | \(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik}\) | 3.9303 | 8 | 4 | 96.07 |
| Conditional on Y | \(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk}\) | 2.9203 | 8 | 4 | 95.06 |
| Conditional on Z | \(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk}\) | 29.729 | 9 | 3 | 123.87 |
| Interaksi | Pengujian | \(\Delta\) deviance | \(\Delta\) df | Chi-square Tabel | Keputusan | Keterangan |
|---|---|---|---|---|---|---|
| XYZ | Saturated vs Homogen | 1.798 | 2 | 5.991 | Tidak Tolak \(H_0\) | tidak ada interaksi |
| YZ | Conditional on X vs Homogen | 2.132 | 2 | 5.991 | Tidak Tolak \(H_0\) | tidak ada interaksi |
| XZ | Conditional on Y vs Homogen | 1.122 | 2 | 5.991 | Tidak Tolak \(H_0\) | tidak ada interaksi |
| XY | Conditional on Z vs Homogen | 27.931 | 1 | 3.841 | Tolak \(H_0\) | ada interaksi |
Berdasarkan hasil di atas, diketahui bahwa asosiasi yang signifikan hanya ditemukan antara tingkat pendidikan (X) dan sikap terhadap teknologi (Y). Sehingga, model terbaik adalah:
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} \]
Model terbaik adalah model log-linear tanpa interaksi tiga arah dan hanya memuat interaksi dua arah antara tingkat pendidikan dan sikap terhadap teknologi.
Model terbaik dipilih berdasarkan hasil pengujian interaksi yang signifikan, yaitu hanya terdapat interaksi dua arah antara tingkat pendidikan (X) dan sikap terhadap teknologi (Y).
Model terbaik dirumuskan sebagai:
\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} \]
# Model Terbaik
bestmodel <- glm(counts ~ x.educ + y.tech + z.job + x.educ*y.tech, family = poisson(link = "log"))
summary(bestmodel)
##
## Call:
## glm(formula = counts ~ x.educ + y.tech + z.job + x.educ * y.tech,
## family = poisson(link = "log"))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.11204 0.08541 48.147 < 2e-16 ***
## x.educRendah -0.32542 0.11510 -2.827 0.004694 **
## y.techSetuju 0.74721 0.09048 8.258 < 2e-16 ***
## z.jobPegawai -0.14108 0.07525 -1.875 0.060805 .
## z.jobProfesional 0.07599 0.07121 1.067 0.285929
## x.educRendah:y.techSetuju 0.44904 0.13493 3.328 0.000875 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 266.348 on 11 degrees of freedom
## Residual deviance: 15.024 on 6 degrees of freedom
## AIC: 101.93
##
## Number of Fisher Scoring iterations: 4
Dari hasil ringkasan model di atas, best model memiliki nilai AIC sebesar 101.93. Meskipun demikian, nilai ini tidak lebih kecil dibandingkan model homogeneous (97.93), conditional on X (96.067), maupun conditional on Y (95.057). Akan tetapi, best model tetap dipilih karena didasarkan pada uji interaksi, di mana hanya interaksi antara tingkat pendidikan (X) dan sikap terhadap teknologi (Y) yang terbukti signifikan, sesuai dengan tujuan pemodelan yang mengutamakan parsimoni dan relevansi hubungan antar variabel.
Untuk memudahkan interpretasi model log-linear, kita mengubah koefisien log menjadi rasio ekspektasi dengan cara mengeksponensialkan nilai koefisien.
# Interpretasi koefisien model terbaik
interpretasi_koef <- data.frame(
koef = bestmodel$coefficients,
exp_koef = exp(bestmodel$coefficients)
)
interpretasi_koef
\(\exp(\lambda^{X}_{\text{Rendah}}) = \exp(-0.32542) \approx 0.722\) → rasio ekspektasi
Tanpa memperhitungkan faktor sikap terhadap teknologi dan pekerjaan, tingkat pendidikan rendah memiliki ekspektasi jumlah sebesar 0.722 kali dibandingkan dengan pendidikan tinggi.
\(\exp(\lambda^{Y}_{\text{Setuju}}) = \exp(0.74721) \approx 2.111\) → rasio ekspektasi
Tanpa memperhitungkan pendidikan dan jenis pekerjaan, responden yang setuju terhadap teknologi memiliki ekspektasi frekuensi sekitar 2.111 kali lebih besar dibandingkan dengan yang tidak setuju.
\(\exp(\lambda^{Z}_{\text{Pegawai}}) = \exp(-0.14108) \approx 0.868\) → rasio ekspektasi
Tanpa memperhitungkan pendidikan dan sikap terhadap teknologi, pekerjaan sebagai pegawai memiliki ekspektasi jumlah 0.868 kali dibandingkan dengan wirausaha.
\(\exp(\lambda^{Z}_{\text{Profesional}}) = \exp(0.07599) \approx 1.079\) → rasio ekspektasi
Tanpa memperhitungkan pendidikan dan sikap terhadap teknologi, pekerjaan sebagai profesional memiliki ekspektasi jumlah sekitar 1.079 kali dibandingkan dengan wirausaha.
\(\exp(\lambda^{XY}_{\text{Rendah,Setuju}}) = \exp(0.44904) \approx 1.567\) → rasio odds
Interaksi menunjukkan bahwa untuk individu dengan pendidikan rendah yang setuju terhadap teknologi, ekspektasi jumlah meningkat sekitar 1.567 kali dibandingkan dengan baseline (pendidikan tinggi dan tidak setuju).
Eksponensial dari koefisien (\(\exp(\beta)\)) dalam model log-linear menunjukkan rasio ekspektasi atau expected count ratio relatif terhadap kategori referensi.
Semua interpretasi diasumsikan ceteris paribus, yaitu faktor lain tetap konstan.
Fitted values dari model log-linear terbaik digunakan untuk melihat seberapa baik model memprediksi jumlah observasi berdasarkan kombinasi tingkat pendidikan, sikap terhadap teknologi, dan jenis pekerjaan.
# Fitted values dari model terbaik
data.frame( Pekerjaan = z.job, Pendidikan = x.educ, Sikap = y.tech, Observasi = counts, Prediksi = bestmodel$fitted.values )
Untuk menghitung fitted value secara manual, dapat menggunakan rumus sebagai berikut:
\[ \hat{\mu}_{ijk} = \exp\left( \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} \right) \]
Berikut adalah perhitungan fitted value secara manual
\[ \begin{aligned} \hat{\mu}_{111} &= \exp(4.265 - 0.593 + 0.483 + 0.01986 + 0.658) \\ &= \exp(4.833) = 125.595 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{112} &= \exp(4.265 - 0.593 + 0.483 + 0.381 + 0.658) \\ &= \exp(5.195) = 180.279 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{113} &= \exp(4.265 - 0.593 + 0.483 + 0 + 0.658) \\ &= \exp(4.813) = 123.126 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{121} &= \exp(4.265 - 0.593 + 0 + 0.01986 + 0) \\ &= \exp(3.692) = 40.109 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{122} &= \exp(4.265 - 0.593 + 0 + 0.381 + 0) \\ &= \exp(4.053) = 57.572 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{123} &= \exp(4.265 - 0.593 + 0 + 0 + 0) \\ &= \exp(3.672) = 39.320 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{211} &= \exp(4.265 + 0 + 0.483 + 0.01986 + 0) \\ &= \exp(4.768) = 117.691 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{212} &= \exp(4.265 + 0 + 0.483 + 0.381 + 0) \\ &= \exp(5.1295) = 168.933 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{213} &= \exp(4.265 + 0 + 0.483 + 0 + 0) \\ &= \exp(4.748) = 115.377 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{221} &= \exp(4.265 + 0 + 0 + 0.01986 + 0) \\ &= \exp(4.285) = 72.605 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{222} &= \exp(4.265 + 0 + 0 + 0.381 + 0) \\ &= \exp(4.646) = 104.217 \end{aligned} \]
\[ \begin{aligned} \hat{\mu}_{223} &= \exp(4.265 + 0 + 0 + 0 + 0) \\ &= \exp(4.265) = 71.178 \end{aligned} \]
Nilai \(\hat{\mu}_{ijk}\) bersifat invariant terhadap perubahan kategori referensi dari variabel prediktor yang digunakan, yang berarti meskipun kategori referensi diganti, hasil prediksi \(\hat{\mu}_{ijk}\) tetap sama. Perubahan referensi hanya mempengaruhi nilai koefisien model, tetapi tidak mengubah nilai fitted value yang dihasilkan.