Puji syukur kami panjatkan kepada Tuhan Yang Maha Esa atas segala rahmat dan karunia-Nya, sehingga buku ini dapat diselesaikan dengan baik. Buku ini disusun untuk memberikan wawasan dan pemahaman yang mendalam tentang berbagai teknik analisis data kategori, yang merupakan alat penting dalam riset dan pengambilan keputusan di berbagai bidang.
Penulisan buku ini berawal dari kebutuhan untuk mengungkapkan secara sistematis dan komprehensif konsep-konsep dasar serta penerapan praktis dari model-model regresi yang sering digunakan dalam penelitian sosial, ekonomi, kesehatan, serta bidang lainnya. Pembaca diharapkan dapat memahami teori dasar dari masing-masing model statistik tersebut, serta aplikasi praktisnya dalam analisis data yang lebih kompleks.
Kami menyadari bahwa buku ini belum sempurna dan masih banyak kekurangan. Oleh karena itu, kritik dan saran yang konstruktif dari pembaca sangat kami harapkan untuk penyempurnaan di masa yang akan datang.
Terima kasih kami ucapkan kepada semua pihak yang telah memberikan dukungan dan bantuan selama penyusunan buku ini. Semoga buku ini dapat memberikan manfaat yang besar bagi pembaca, baik yang baru memulai belajar statistik maupun yang sudah berpengalaman, dalam memahami dan mengaplikasikan teknik analisis regresi dalam penelitian mereka.
Akhir kata, kami berharap buku ini dapat menjadi referensi yang bermanfaat dan memberikan kontribusi bagi perkembangan ilmu pengetahuan dan praktik analisis statistik.
Penulis
Dekhsa Khoiriyatul Fajri
Sumedang, 11 Juni 2025
Analisis data kategori adalah proses pengolahan dan interpretasi data yang terdiri dari kategori atau label untuk mengidentifikasi pola, hubungan, atau tren yang signifikan. Tujuan dari analisis ini adalah untuk menghasilkan wawasan yang dapat digunakan untuk pengambilan keputusan, pengembangan kebijakan, atau pengujian hipotesis dalam berbagai konteks. Dengan menggunakan teknik seperti tabulasi silang, uji statistik, atau model prediktif, analisis data kategori membantu mengubah data kualitatif menjadi informasi yang terukur dan bermakna, sehingga mendukung pemahaman yang lebih baik terhadap fenomena yang sedang dipelajari.
Buku ini hadir untuk mengisi tujuan tersebut dengan menyediakan panduan praktis dan terstruktur tentang analisis data kategori. Tujuan utama buku ini adalah untuk membekali pembaca dengan pemahaman konseptual yang kuat serta keterampilan teknis untuk mengolah dan menganalisis data kategori menggunakan pendekatan yang relevan dan terkini. Buku ini dirancang untuk menjadi sumber referensi baik bagi pemula yang baru mengenal analisis data kategori maupun bagi profesional yang ingin memperdalam keahlian mereka dalam metode analisis khusus.
Cakupan buku ini meliputi pengenalan konsep dasar data kategori, teknik visualisasi data, metode statistik untuk analisis data kategori seperti uji chi-square, analisis regresi logistik, analisis Log Linear, dan lain sebagainya menggunakan R. Selain itu, buku ini juga membahas aplikasi praktis analisis data kategori dalam berbagai konteks, termasuk studi kasus nyata yang relevan. Dengan pendekatan yang berorientasi pada praktik dan contoh-contoh yang jelas, buku ini diharapkan dapat menjadi panduan yang bermanfaat bagi siapa saja yang ingin memahami dan memanfaatkan potensi data kategori secara efektif.
Analisis data kategori merupakan cabang statistik yang digunakan untuk menganalisis data yang bersifat kategorik, yaitu data yang dapat dikelompokkan ke dalam kategori atau kelas tertentu, baik nominal (tanpa urutan) maupun ordinal (dengan urutan). Data kategorik sering kali muncul dalam berbagai bidang seperti sosiologi, kesehatan, dan pemasaran, misalnya dalam bentuk status kelulusan (lulus/gagal), jenis kelamin (pria/wanita), atau tingkat kepuasan (sangat puas, puas, tidak puas). Tujuan utama analisis ini adalah untuk memahami hubungan antar variabel kategorik, memodelkan probabilitas kejadian, dan membuat inferensi statistik.
Nominal: Tidak memiliki urutan (misalnya, jenis kelamin).
Ordinal: Memiliki urutan (misalnya, tingkat pendidikan).
Biner: Dua kategori (misalnya, ya/tidak).
Multikategori: Lebih dari dua kategori.
Data kategori tidak diukur dalam skala numerik kontinu dan biasanya dianalisis dengan frekuensi dan proporsi.
Segmentasi pelanggan berdasarkan kategori seperti usia, jenis kelamin, atau preferensi produk.
Analisis loyalitas pelanggan: apakah pelanggan dengan kategori tertentu (misalnya pelanggan lama vs baru) lebih cenderung melakukan pembelian ulang.
Menentukan strategi promosi yang tepat berdasarkan kategori demografi.
Mengidentifikasi faktor risiko penyakit berdasarkan kategori seperti jenis kelamin, gaya hidup (merokok/tidak), dan pekerjaan.
Analisis efektivitas pengobatan berdasarkan kelompok pasien (misalnya berdasarkan usia atau riwayat medis).
Mengklasifikasikan jenis penyakit dan melihat penyebarannya di kelompok populasi tertentu.
Menganalisis hubungan antara latar belakang siswa (jenis kelamin, asal sekolah, jurusan) dengan prestasi akademik.
Menentukan strategi pengajaran yang efektif berdasarkan kategori gaya belajar siswa.
Evaluasi program pendidikan berdasarkan kategori tingkat kepuasan siswa.
Menilai kesejahteraan masyarakat berdasarkan kategori pekerjaan, status pernikahan, atau tempat tinggal.
Menganalisis respon publik terhadap kebijakan pemerintah berdasarkan kategori usia atau wilayah.
Klasifikasi kelompok rentan untuk penyaluran bantuan sosial.
Menguji hubungan antar variabel kategorik (contoh: apakah ada hubungan antara jenis kelamin dan preferensi media sosial).
Digunakan dalam analisis statistik inferensial seperti uji chi-square.
Membantu membangun model prediktif dengan data input kategorik.
Digunakan untuk eksperimen biner, yaitu suatu percobaan yang hanya memiliki dua hasil (sukses atau gagal). Fungsi probabilitasnya:
\[ P(X = x) = p^x (1 - p)^{1 - x},\quad x \in \{0, 1\} \]
Keterangan Notasi:
𝑋: Variabel acak biner (0 atau 1)
𝑝: Probabilitas sukses (X = 1)
Contoh: Seorang peneliti ingin mensimulasikan data hasil ujian dari 20 siswa, di mana setiap siswa hanya memiliki dua kemungkinan: 1 jika lulus dan 0 jika tidak lulus
Peluang kelulusan setiap siswa adalah 0.7. Hitung banyaknya siswa yang lulus!
Gunakan distribusi binomial dengan ukuran 1 (size = 1)
untuk mensimulasikan distribusi Bernoulli.
library(knitr)
## Warning: package 'knitr' was built under R version 4.3.3
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.3.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
##
## group_rows
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
set.seed(123)
n <- 20
p <- 0.7
hasil_ujian <- rbinom(n = n, size = 1, prob = p)
hasil_ujian
## [1] 1 0 1 0 0 1 1 0 1 1 0 1 1 1 1 0 1 1 1 0
Kesimpulan: berdasarkan perhitungan di atas, dapat diketahui bahwa jumlah siswa yang lulus (1) sebanyak 13 siswa.
Distribusi binomial adalah jumlah keberhasilan dalam n percobaan Bernoulli.
\[ P(X = k) = \binom{n}{k} p^k (1 - p)^{n - k} \] \(P(X = k)\): Probabilitas bahwa terdapat tepat \(k\) keberhasilan dari \(n\) percobaan.
\(n\): Jumlah percobaan (banyaknya uji atau sampel yang diamati).
\(k\): Jumlah keberhasilan yang diinginkan.
\(p\): Probabilitas keberhasilan dalam satu percobaan.
\(1 - p\): Probabilitas kegagalan dalam satu percobaan.
\(\binom{n}{k}\): Kombinasi, yaitu banyak cara memilih \(k\) keberhasilan dari \(n\) percobaan, yang dirumuskan sebagai:
\[ \binom{n}{k} = \frac{n!}{k!(n - k)!} \] Contoh:
Sebuah perusahaan ingin mengetahui seberapa sering produknya berhasil lulus uji kualitas. Probabilitas sebuah produk lulus uji adalah 0.8. Jika diuji 10 produk secara acak, berapa kemungkinan ada tepat 8 produk yang lulus? Simulasikan 10 kali percobaan pengujian 10 produk menggunakan distribusi binomial, dan hitung berapa kali muncul tepat 8 produk yang lulus.
set.seed(123)
# Parameter binomial
n <- 10 # Jumlah produk diuji
p <- 0.8 # Probabilitas lulus uji
simulasi <- 10 # Jumlah simulasi
# Simulasi distribusi binomial
hasil <- rbinom(simulasi, size = n, prob = p)
# Menghitung berapa kali muncul tepat 8 produk yang lulus
jumlah_8_lulus <- sum(hasil == 8)
jumlah_8_lulus
## [1] 4
Kesimpulan: Berdasarkan perhitungan di atas, dapat diketahui bahwa dari 10 kali percobaan, sebanyak 3 kali terjadi kejadian tepat 8 produk lulus uji.
Distribusi multinomial adalah generalisasi dari distribusi binomial, yang digunakan untuk menggambarkan data kategorik dengan lebih dari dua kategori. Misalkan terdapat \(k\) kategori, dan setiap observasi hanya dapat berada di satu kategori. Jika kita melakukan \(n\) percobaan independen, dengan probabilitas masing-masing kategori \(\pi_1, \pi_2, \ldots, \pi_k\) (di mana \(\sum_{j=1}^k \pi_j = 1\)), maka distribusi multinomial didefinisikan sebagai: \[P(Y_1 = y_1, Y_2 = y_2, \ldots, Y_k = y_k) = \frac{n!}{y_1! y_2! \cdots y_k!} \pi_1^{y_1} \pi_2^{y_2} \cdots \pi_k^{y_k}\] Di mana:
\(Y_j\): Jumlah observasi dalam kategori ke-\(j\).
\(y_j\): Nilai observasi untuk kategori ke-\(j\), dengan \(\sum_{j=1}^k y_j = n\).
\(\pi_j\): Probabilitas kategori ke-\(j\).
Contoh soal: Sebuah survei dilakukan terhadap 100 orang mengenai jenis minuman favorit mereka, dengan tiga pilihan: Teh, Kopi, dan Jus. Diketahui bahwa probabilitas seseorang memilih Teh adalah 0.4, Kopi 0.35, dan Jus 0.25. Simulasikan hasil survei tersebut menggunakan distribusi multinomial sebanyak 100 responden, dan jawab pertanyaan berikut:
Jawaban:
set.seed(123)
# Parameter simulasi
n <- 100 # jumlah responden
prob <- c(Teh = 0.4, Kopi = 0.35, Jus = 0.25)
# Simulasi multinomial
hasil <- rmultinom(n = 1, size = n, prob = prob)
hasil
## [,1]
## Teh 39
## Kopi 36
## Jus 25
Kesimpulan: Berdasarkan hasil perhitungan, dapat diketahui bahwa sekitar 39 orang memilih teh, sekitar 36 orang memilih kopi, sekitar 25 orang memilih jus.
Distribusi poisson digunakan untuk menghitung jumlah kejadian dalam interval tertentu. Fungsi probabilitasnya:
\[ P(X = k) = \frac{e^{-\lambda} \lambda^k}{k!} \]
𝑋: jumlah kejadian dalam interval tertentu
𝜆: rata-rata kejadian dalam interval tersebut
𝑘: jumlah kejadian yang diamati
Contoh soal: ## Soal
Sebuah call center menerima rata-rata 4 telepon per jam. Misalkan jumlah telepon yang masuk mengikuti distribusi Poisson. Buatlah simulasi sebanyak 10 data untuk merepresentasikan jumlah telepon per jam yang masuk ke call center tersebut.
Jawaban
# Set seed untuk replikasi
set.seed(123)
# Simulasi 100 data dari distribusi Poisson dengan lambda = 4
telepon <- rpois(10, lambda = 4)
telepon
## [1] 3 6 3 6 7 1 4 7 4 4
Kesimpulan: Berdasarkan perhitungan di atas, dapat diketahui bahwa dalam 10 jam, jumlah kejadian panggilan telepon masuk yang terjadi pada satu jam pertama adalah sebanyak 4, pada dua jam kedua adalah sebanyak 3, pada tiga jam pertama berjumlah 4, dan seterusnya.
Prospective sampling melibatkan pengumpulan data secara prospektif, di mana peneliti mengamati subjek dari waktu saat ini ke masa depan untuk mencatat kejadian atau hasil tertentu. Pendekatan ini sering digunakan dalam studi yang dirancang untuk mengevaluasi hubungan sebab-akibat.
Desain eksperimen melibatkan intervensi aktif oleh peneliti, di mana subjek secara acak ditugaskan ke kelompok perlakuan atau kontrol untuk mengamati efek dari intervensi tersebut. Dalam analisis data kategorik, eksperimen sering digunakan untuk mengevaluasi efektivitas suatu perlakuan terhadap hasil kategorik, seperti keberhasilan pengobatan (sembuh/tidak sembuh).
Karakteristik:
Contoh: Sebuah uji klinis dilakukan untuk menguji efektivitas obat baru terhadap penyakit tertentu. Pasien secara acak diberi obat atau plasebo, dan setelah periode tertentu, status kesehatan mereka (sembuh/tidak sembuh) dicatat.
Studi kohort melibatkan pengamatan sekelompok individu yang memiliki karakteristik tertentu (kohort) selama periode waktu tertentu untuk mencatat kejadian atau hasil tertentu. Dalam data kategorik, studi ini sering digunakan untuk menganalisis hubungan antara paparan tertentu dan hasil kategorik.
Karakteristik:
Contoh: Sebuah studi kohort dilakukan untuk meneliti hubungan antara kebiasaan merokok dan kejadian kanker paru-paru. Sejumlah individu (perokok dan bukan perokok) dipantau selama 10 tahun, dan kejadian kanker paru-paru (ada/tidak ada) dicatat.
Retrospective sampling melibatkan pengumpulan data secara retrospektif, di mana peneliti melihat ke masa lalu untuk mengidentifikasi subjek berdasarkan hasil atau status tertentu, kemudian menelusuri faktor risiko atau paparan yang mungkin berkontribusi.
Studi kasus-kontrol dimulai dengan mengidentifikasi individu yang memiliki hasil tertentu (kasus) dan individu yang tidak memiliki hasil tersebut (kontrol), kemudian menelusuri paparan masa lalu untuk mengevaluasi hubungannya dengan hasil tersebut. Studi ini sering digunakan dalam analisis data kategorik untuk meneliti penyebab langka.
Karakteristik:
Contoh: Sebuah studi kasus-kontrol dilakukan untuk meneliti hubungan antara paparan asbes dan mesothelioma. Peneliti mengambil sampel pasien mesothelioma (kasus) dan individu sehat (kontrol), lalu menelusuri riwayat paparan asbes mereka.
Studi kohort retrospektif mirip dengan studi kohort, tetapi data dikumpulkan dari catatan masa lalu. Peneliti memilih kohort berdasarkan status paparan di masa lalu, kemudian menelusuri hasil yang telah terjadi hingga saat ini.
Karakteristik:
Contoh: Sebuah studi kohort retrospektif dilakukan untuk meneliti hubungan antara paparan radiasi di tempat kerja dan kejadian leukemia. Peneliti menggunakan catatan karyawan dari 20 tahun lalu untuk mengidentifikasi kelompok yang terpapar radiasi, lalu memeriksa catatan medis untuk mencatat kejadian leukemia.
Dalam memilih desain sampling untuk analisis data kategorik, beberapa hal perlu diperhatikan:
Arah Waktu: Prospective sampling cocok untuk mengevaluasi hubungan sebab-akibat secara langsung, sedangkan retrospective sampling lebih efisien untuk meneliti kejadian langka atau ketika data historis tersedia.
Bias dan Validitas: Retrospective sampling rentan terhadap bias recall (misalnya, subjek mungkin salah mengingat paparan masa lalu). Prospective sampling dapat mengurangi bias ini, tetapi memerlukan waktu dan biaya lebih besar.
Efisiensi: Studi kasus-kontrol lebih efisien untuk kejadian langka, sedangkan studi kohort lebih cocok untuk paparan yang umum.
Tabel kontngensi adalah tabel klasifikasi silang dari variabel-variabel yang bersifat kategori. Tabel kontingensi dapat menggambarkan hubungan diantara variabel-variabel kategori tersebut. Terdapat beberapa jenis tabel kontingensi, salah satunya adalah tabel kontingensi dua arah. Tabel kontingensi dua arah adalah tabel kontingensi yang memiliki dua variabel kategori yang terlibat. Sedangkan, tabel kontingensi 2x2 adalah tabel kontingensi yang memiliki dua variabel yang mana setiap variabel tersebut memiliki 2 kategori.
| Kategori 1 (+) | Kategori 2 (-) | Total | |
|---|---|---|---|
| Grup 1 | \(n_{11}\) | \(n_{12}\) | \(n_{1.}\) |
| Grup 2 | \(n_{21}\) | \(n_{22}\) | \(n_{2.}\) |
| Total | \(n_{.1}\) | \(n_{.2}\) | \(n\) |
Distribusi peluang dari tabel kontingensi terdiri dari:
Distribusi peluang bersama (Joint probability)
\[ P(X=i, Y=j) = \frac{n_{ij}}{n} \]
Distribusi peluang marginal (Marginal probability)
\[ P(X=i) = \sum_{j} P(X=i, Y=j) \]
\[ P(Y=j) = \sum_{i} P(X=i, Y=j) \]
Distribusi peluang bersyarat (Conditional probability)
\[ P(Y=j | X=i) = \frac{P(X=i, Y=j)}{P(X=i)} \]
\[ P(X=i | Y=j) = \frac{P(X=i, Y=j)}{P(Y=j)} \]
Ukuran Asosiasi adalah statistik yang menyediakan standar untuk menilai hubungan antara variabel yang diamati dalam tabel kontingensi. Terdapat tiga jenis ukuran asosiasi yang dapat digunakan untuk tabel kontingensi dua arah yaitu,
Beda Peluang (Difference Probability)
Beda peluang adalah selisih dari variabel respons dengan syarat variabel prediktor yang berbeda. Beda peluang memberikan indikasi ada tidaknya asosiasi antara keduanya.
\[ RD = P(Y=1|X=1) - P(Y=1|X=0) \]
Resiko Relatif (Relative Risk)
Resiko relatif (RR) adalah perbandingan antara risiko subjek yang terpapar dan tidak terpapar. Resiko relatif biasanya digunakan dalam studi prospektif seperti studi kohort.
\[ RR = \frac{P(Y=1|X=1)}{P(Y=1|X=0)} \]
Rasio Odds (Odds Ratio)
Odds ratio adalah perbandingan dari dua variabel respons terhadap variabel prediktor yang sama. Odds ratio biasanya digunakan dalam studi kasus kontrol.
\[ OR = \frac{P(Y=1|X=1) / P(Y=0|X=1)}{P(Y=1|X=0) / P(Y=0|X=0)} \]
Contoh Kasus Tabel Kontingensi Dua Arah:
Misalkan kita memiliki data hasil survei mengenai kebiasaan merokok dan status kesehatan jantung. Hitung peluang bersama, peluang marjinal, peluang bersyarat, dan ukuran asosiasinya!
| Sakit Jantung (Y=1) | Tidak Sakit Jantung (Y=0) | Total | |
|---|---|---|---|
| Merokok (X=1) | 30 | 10 | 40 |
| Tidak Merokok (X=0) | 20 | 40 | 60 |
| Total | 50 | 50 | 100 |
Peluang Bersama
Peluang bersama adalah peluang kejadian gabungan dari dua variabel, dihitung sebagai:
\(P(X=1, Y=1)\) (Merokok dan Sakit Jantung):
\[ P(X=1, Y=1) = \frac{n_{11}}{n} = \frac{30}{100} = 0.3 \]
\(P(X=1, Y=0)\) (Merokok dan Tidak Sakit Jantung):
\[ P(X=1, Y=0) = \frac{n_{12}}{n} = \frac{10}{100} = 0.1 \]
\(P(X=0, Y=1)\) (Tidak Merokok dan Sakit Jantung):
\[ P(X=0, Y=1) = \frac{n_{21}}{n} = \frac{20}{100} = 0.2 \]
\(P(X=0, Y=0)\) (Tidak Merokok dan Tidak Sakit Jantung):
\[ P(X=0, Y=0) = \frac{n_{22}}{n} = \frac{40}{100} = 0.4 \]
Peluang Marjinal (Marginal Probability)
\(P(X=1)\) (Merokok):
\[ P(X=1) = \frac{n_{1.}}{n} = \frac{40}{100} = 0.4 \]
\(P(X=0)\) (Tidak Merokok):
\[ P(X=0) = \frac{n_{2.}}{n} = \frac{60}{100} = 0.6 \]
\(P(Y=1)\) (Sakit Jantung):
\[ P(Y=1) = \frac{n_{.1}}{n} = \frac{50}{100} = 0.5 \]
\(P(Y=0)\) (Tidak Sakit Jantung):
\[ P(Y=0) = \frac{n_{.2}}{n} = \frac{50}{100} = 0.5 \]
Peluang Bersyarat (Conditional Probability)
\(P(Y=1 | X=1)\) (Sakit Jantung diberikan Merokok):
\[ P(Y=1 | X=1) = \frac{n_{11}}{n_{1.}} = \frac{30}{40} = 0.75 \]
\(P(Y=1 | X=0)\) (Sakit Jantung diberikan Tidak Merokok):
\[ P(Y=1 | X=0) = \frac{n_{21}}{n_{2.}} = \frac{20}{60} = 0.3333 \]
\(P(X=1 | Y=1)\) (Merokok diberikan Sakit Jantung):
\[ P(X=1 | Y=1) = \frac{n_{11}}{n_{.1}} = \frac{30}{50} = 0.6 \]
\(P(X=1 | Y=0)\) (Merokok diberikan Tidak Sakit Jantung):
\[ P(X=1 | Y=0) = \frac{n_{12}}{n_{.2}} = \frac{10}{50} = 0.2 \]
Perhitungan R
Peluang Bersama
n <- 100
P_X1_Y1 <- 30 / n
P_X1_Y0 <- 10 / n
P_X0_Y1 <- 20 / n
P_X0_Y0 <- 40 / n
P_X1_Y1; P_X1_Y0; P_X0_Y1; P_X0_Y0
## [1] 0.3
## [1] 0.1
## [1] 0.2
## [1] 0.4
Peluang Marjinal
P_X1 <- (30 + 10) / n
P_X0 <- (20 + 40) / n
P_Y1 <- (30 + 20) / n
P_Y0 <- (10 + 40) / n
P_X1; P_X0; P_Y1; P_Y0
## [1] 0.4
## [1] 0.6
## [1] 0.5
## [1] 0.5
Peluang Bersyarat
P_Y1_given_X1 <- P_X1_Y1 / P_X1
P_Y0_given_X1 <- P_X1_Y0 / P_X1
P_Y1_given_X0 <- P_X0_Y1 / P_X0
P_Y0_given_X0 <- P_X0_Y0 / P_X0
P_X1_given_Y1 <- P_X1_Y1 / P_Y1
P_X0_given_Y1 <- P_X0_Y1 / P_Y1
P_X1_given_Y0 <- P_X1_Y0 / P_Y0
P_X0_given_Y0 <- P_X0_Y0 / P_Y0
P_Y1_given_X1; P_Y0_given_X1; P_Y1_given_X0; P_Y0_given_X0
## [1] 0.75
## [1] 0.25
## [1] 0.3333333
## [1] 0.6666667
P_X1_given_Y1; P_X0_given_Y1; P_X1_given_Y0; P_X0_given_Y0
## [1] 0.6
## [1] 0.4
## [1] 0.2
## [1] 0.8
Kesimpulan
Ukuran Asosiasi
Di mana \(\hat{p}{1|1} = P(Y=1 | X=1) = 0.75\) dan \(\hat{p}{1|2} = P(Y=1 | X=0) = 0.3333\).
\[ \hat{\text{RD}} = 0.75 - 0.3333 = 0.4167 \]
Relative Risk (RR)
\[ \hat{\text{RR}} = \frac{0.75}{0.3333} = 2.25 \]
Odds Ratio (OR)
\[ \hat{\text{OR}} = \frac{30 \times 40}{10 \times 20} = \frac{1200}{200} = 6 \]
Perhitungan R
Beda Peluang
RD <- P_Y1_given_X1 - P_Y1_given_X0
RD
## [1] 0.4166667
Resiko Relatif (Relative Risk)
RR <- P_Y1_given_X1 / P_Y1_given_X0
RR
## [1] 2.25
Odds Ratio (OR)
OR <- (P_Y1_given_X1 / P_Y0_given_X1) / (P_Y1_given_X0 / P_Y0_given_X0)
OR
## [1] 6
Kesimpulan:
Kesimpulan: Karena nilai RR dan OR jauh lebih besar dari 1, maka merokok dapat dianggap sebagai faktor risiko utama terhadap penyakit jantung.
Estimasi titik bertujuan untuk memperkirakan satu nilai spesifik dari suatu parameter. Estimasi titik proporsi dihitung dengan rumus:
\[ \hat{p} = \frac{x}{n} \]
dimana:
Estimasi interval digunakan untuk memberikan rentang nilai yang diyakini memiliki nilai parameter populasi dengan tingkat kepercayaan tertentu (biasanya 99% atau 95%). Rumus untuk menghitung interval kepercayaan untuk proporsi adalah:
\[ \hat{p} \pm Z_{\alpha/2} \sqrt{ \frac{ \hat{p}(1 - \hat{p}) }{n} } \]
dimana:
Uji proporsi digunakan untuk menguji perbedaan proporsi antara dua kelompok dalam tabel kontingensi 2×2.
Rumus Uji Proporsi Dua Sampel:
\[ Z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{\hat{p}(1 - \hat{p})\left(\frac{1}{n_1} + \frac{1}{n_2}\right)}} \] di mana:
\(\hat{p}_1 = \frac{n_{11}}{n_1}\) (proporsi kelompok 1)
\(\hat{p}_2 = \frac{n_{21}}{n_2}\) (proporsi kelompok 2)
\(\hat{p} = \frac{n_{11} + n_{21}}{n_1 + n_2}\) (proporsi gabungan).
Risk different digunakan untuk mengukur selisih absolut risiko antara
dua kelompok.
\[
RD = \hat{p}_1 - \hat{p}_2
\]
Standar Error:
\[
SE(RD) = \sqrt{\frac{\hat{p}_1(1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2(1
- \hat{p}_2)}{n_2}}
\] Statistik Uji \(z\):
\[ z = \frac{\hat{\text{RD}}}{\hat{\text{SE RD}}} \]
Relative Risk digunakan untuk melakukan perbandingan risiko antara
kelompok terpapar dan tidak terpapar.
\[
RR = \frac{\hat{p}_1}{\hat{p}_2}
\]
Standar Error log(RR):
\[
SE(\ln RR) = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_1} + \frac{1}{n_{21}} +
\frac{1}{n_2}}
\]
Uji statistik z:\[ z = \frac{\log(\hat{\text{RR}})}{\hat{\text{SE log RR}}} \]
Odds ratio digunakan untuk melakukan perbandingan odds antara
kelompok terpapar dan tidak terpapar.
\[
OR = \frac{n_{11} \times n_{22}}{n_{12} \times n_{21}}
\]
Standar Error log(OR):
\[
SE(\ln OR) = \sqrt{\frac{1}{n_{11}} - \frac{1}{n_{12}} +
\frac{1}{n_{21}} - \frac{1}{n_{22}}}
\] Uji statistik z: \[
z = \frac{\log(\hat{\text{OR}})}{\hat{\text{SE log OR}}}
\]
Uji Chi-Square digunakan untuk menguji hubungan statistik antara dua
variabel kategori.
\[
\chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}}
\]
di mana \(E_{ij} = \frac{(n_{i+} \times
n_{+j})}{n}\).
Contoh soal:
Sebuah survei dilakukan untuk mengetahui apakah ada hubungan antara jenis kelamin (Pria/Wanita) dan preferensi terhadap produk A (Suka/Tidak Suka).
Berikut adalah hasil data dalam bentuk tabel kontingensi 2x2:
| Suka | Tidak Suka | Total | |
|---|---|---|---|
| Pria | 50 | 10 | 60 |
| Wanita | 20 | 40 | 60 |
| Total | 70 | 50 | 120 |
Tentukan:
Apakah jenis kelamin dan preferensi produk A saling bebas (independen)?
Apakah proporsi orang yang menyukai produk A berbeda secara signifikan antara pria dan wanita?
Apakah terdapat asosiasi antara jenis kelamin dan preferensi produk?
Hipotesis:
Tingkat Signifikansi:
alfa: 5% = 0,05
Perhitungan manual
Hitung frekuensi yang diharapkan (\(e_{ij}\)):
\(e_{11}\) (Pria, Suka):
\[ e_{11} = \frac{n_{1.} \cdot n_{.1}}{n} = \frac{60 \times 70}{120} = \frac{4200}{120} = 35 \]
\(e_{12}\) (Pria, Tidak Suka):
\[ e_{12} = \frac{n_{1.} \cdot n_{.2}}{n} = \frac{60 \times 50}{120} = \frac{3000}{120} = 25 \]
\(e_{21}\) (Wanita, Suka):
\[ e_{21} = \frac{n_{2.} \cdot n_{.1}}{n} = \frac{60 \times 70}{120} = \frac{4200}{120} = 35 \]
\(e_{22}\) (Wanita, Tidak Suka):
\[ e_{22} = \frac{n_{2.} \cdot n_{.2}}{n} = \frac{60 \times 50}{120} = \frac{3000}{120} = 25 \]
Hitung statistik Chi-Square:
\[ \chi^2 = \frac{(n_{11} - e_{11})^2}{e_{11}} + \frac{(n_{12} - e_{12})^2}{e_{12}} + \frac{(n_{21} - e_{21})^2}{e_{21}} + \frac{(n_{22} - e_{22})^2}{e_{22}} \]
\[ \chi^2 = \frac{(50 - 35)^2}{35} + \frac{(10 - 25)^2}{25} + \frac{(20 - 35)^2}{35} + \frac{(40 - 25)^2}{25} \]
\[ \chi^2 = \frac{(15)^2}{35} + \frac{(-15)^2}{25} + \frac{(-15)^2}{35} + \frac{(15)^2}{25} \]
\[ \chi^2 = \frac{225}{35} + \frac{225}{25} + \frac{225}{35} + \frac{225}{25} \]
\[ \chi^2 = \frac{225}{35} + 9 + \frac{225}{35} + 9 \]
\[ \frac{225}{35} = 6.4286 \]
\[ \chi^2 = 6.4286 + 9 + 6.4286 + 9 = 30.8572 \]
Perhitungan Manual di R:
# Data Observasi
observed <- matrix(c(50, 10,
20, 40),
nrow = 2,
byrow = TRUE)
colnames(observed) <- c("Suka", "Tidak Suka")
rownames(observed) <- c("Pria", "Wanita")
observed
## Suka Tidak Suka
## Pria 50 10
## Wanita 20 40
# Total baris dan kolom
row_totals <- rowSums(observed)
col_totals <- colSums(observed)
total <- sum(observed)
# Frekuensi harapan
expected <- outer(row_totals, col_totals) / total
expected
## Suka Tidak Suka
## Pria 35 25
## Wanita 35 25
# Chi-square manual
chi_square_manual <- sum((observed - expected)^2 / expected)
chi_square_manual
## [1] 30.85714
Perhitungan Dengan Package:
# Menggunakan chisq.test
chi_test <- chisq.test(observed, correct = FALSE)
chi_test
##
## Pearson's Chi-squared test
##
## data: observed
## X-squared = 30.857, df = 1, p-value = 2.777e-08
Kriteria Uji:
Tolak H0 jika nilai χ² > nilai χ² tabel atau p-value < alfa
Keputusan:
Karena χ² (30,86) > χ² tabel (3,84) dan p-value (2,777e-08) < alfa (0,05), maka tolak H0. Artinya, terdapat hubungan antara variabel jenis kelamin terhadap preferensi kesukaan produk A.
Hipotesis
H0: \(p_1 =
p_2\)
(Tidak terdapat perbedaan proporsi antara pria dan wanita yang menyukai
produk A)
H1: \(p_1 \ne
p_2\)
(Terdapat perbedaan proporsi antara pria dan wanita yang menyukai produk
A)
Taraf Signifikansi
alfa = 5% = 0,05
Perhitungan manual
Hitung proporsi:
\[ \hat{p}{1|1} = \frac{n{11}}{n*{1.}} = \frac{50}{60} = 0.8333 \]
\[ \hat{p}{1|2} = \frac{n{21}}{n*{2.}} = \frac{20}{60} = 0.3333 \]
Hitung proporsi gabungan:
\[ \hat{p} = \frac{n*{11} + n*{21}}{n} = \frac{50 + 20}{120} = \frac{70}{120} = 0.5833 \]
Hitung standar error:
\[ \hat{p}(1 - \hat{p}) = 0.5833 \times (1 - 0.5833) = 0.5833 \times 0.4167 = 0.2431 \]
\[ \frac{1}{n_{1.}} + \frac{1}{n_{2.}} = \frac{1}{60} + \frac{1}{60} = \frac{2}{60} = 0.0333 \]
\[ \hat{p}(1 - \hat{p}) \left( \frac{1}{n_{1.}} + \frac{1}{n_{2.}} \right) = 0.2431 \times 0.0333 = 0.0081 \]
\[ \sqrt{0.0081} = 0.09 \]
Hitung statistik \(z\):
\[ z = \frac{\hat{p}{1|1} - \hat{p}{1|2}}{\sqrt{\hat{p}(1 - \hat{p}) \left( \frac{1}{n*{1.}} + \frac{1}{n*{2.}} \right)}} = \frac{0.8333 - 0.3333}{0.09} = \frac{0.5}{0.09} = 5.5556 \]
Perhitungan Manual R
# Data
x1 <- 50 # pria yang suka
n1 <- 60 # total pria
x2 <- 20 # wanita yang suka
n2 <- 60 # total wanita
# Proporsi
p1 <- x1 / n1
p2 <- x2 / n2
p_comb <- (x1 + x2) / (n1 + n2)
# Statistik uji Z
z <- (p1 - p2) / sqrt(p_comb * (1 - p_comb) * (1/n1 + 1/n2))
z
## [1] 5.554921
Perhitungan Dengan Package
prop.test(x = c(50, 20), n = c(60, 60), correct = FALSE)
##
## 2-sample test for equality of proportions without continuity correction
##
## data: c(50, 20) out of c(60, 60)
## X-squared = 30.857, df = 1, p-value = 2.777e-08
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.3479477 0.6520523
## sample estimates:
## prop 1 prop 2
## 0.8333333 0.3333333
Kriteria Uji
Tolak H0 jika nilai Z > nilai Z tabel atau p-value < alfa
Keputusan
Karena Z (5,55) > Z tabel (1,96) dan p-value (2,777e-08) < alfa (0,05), maka tolak H0. Artinya, terdapat perbedaan proporsi antara pria dan wanita yang menyukai produk A.
H0: Tidak ada asosiasi antara jenis kelamin dan preferensi terhadap produk A
H1: Ada asosiasi antara jenis kelamin dan preferensi terhadap produk A
Tingkat Signifikansi
alfa = 5% = 0,05
Perhitungan manual
Hitung \(\log(\hat{\text{OR}})\):
\[ \hat{\text{OR}} = 10 \]
\[ \log(10) \approx 2.3026 \quad (\text{menggunakan logaritma natural}) \]
Hitung standar error:
\[ \frac{1}{n*{11}} + \frac{1}{n*{12}} + \frac{1}{n*{21}} + \frac{1}{n*{22}} = \frac{1}{50} + \frac{1}{10} + \frac{1}{20} + \frac{1}{40} \]
\[ \frac{1}{50} = 0.02, \quad \frac{1}{10} = 0.1, \quad \frac{1}{20} = 0.05, \quad \frac{1}{40} = 0.025 \]
\[ 0.02 + 0.1 + 0.05 + 0.025 = 0.195 \]
\[ \sqrt{0.195} \approx 0.4416 \]
Hitung statistik \(z\):
\[ z = \frac{\log(\hat{\text{OR}})}{\sqrt{\frac{1}{n*{11}} + \frac{1}{n*{12}} + \frac{1}{n*{21}} + \frac{1}{n*{22}}}} = \frac{2.3026}{0.4416} \approx 5.2156 \]
Perhitungan R
# Data
a <- 50 # pria suka
b <- 10 # pria tidak suka
c <- 20 # wanita suka
d <- 40 # wanita tidak suka
n1 <- a + b # total pria
n2 <- c + d # total wanita
p1 <- a / n1
p2 <- c / n2
# Risk Difference (BD)
bd <- p1 - p2
# Standard Error BD
se_bd <- sqrt((p1 * (1 - p1)) / n1 + (p2 * (1 - p2)) / n2)
# Z BD
z_bd <- bd / se_bd
# Hasil
bd; se_bd; z_bd
## [1] 0.5
## [1] 0.07757911
## [1] 6.445034
Relative Risk (RR)
rr <- p1 / p2
se_log_rr <- sqrt((1/a) - (1/n1) + (1/c) - (1/n2))
z_rr <- log(rr) / se_log_rr
# Hasil
rr; se_log_rr; z_rr
## [1] 2.5
## [1] 0.1914854
## [1] 4.785172
Odds Ratio (OR)
or <- (a * d) / (b * c)
se_log_or <- sqrt(1/a + 1/b + 1/c + 1/d)
z_or <- log(or) / se_log_or
# Hasil
or; se_log_or; z_or
## [1] 10
## [1] 0.441588
## [1] 5.214328
Kriteria Uji
Tolak H0 jika nilai Z > nilai Z tabel
Keputusan
Karena nilai dari Z BD (6,44), Z RR (6,79), dan Z OR (5,21) > Z tabel (1,96), maka tolak H0. Artinya, terdapat asosiasi antara jenis kelamin dan preferensi terhadap produk A.
Kesimpulan
Berdasarkan hasil uji Chi-Square, dapat disimpulkan bahwa jenis kelamin berhubungan dengan preferensi terhadap produk A. Dengan kata lain, pria dan wanita punya kecenderungan yang berbeda dalam menilai produk A. Jadi, jenis kelamin bisa jadi faktor yang memengaruhi suka atau tidaknya seseorang terhadap produk tersebut.
Berdasarkan uji proporsi, dapat disimpulkann bahwa terdapat perbedaan nyata antara pria dan wanita dalam menyukai produk A. Pria cenderung lebih menyukai produk A dibandingkan dengan wanita. Jadi, jenis kelamin berpengaruh terhadap preferensi produk A.
Uji asosiasi beda peluang (RD): Berdasarkan hasil perhitungan, nilai RD menunjukkan bahwa pria menyukai produk A jauh lebih tinggi dibanding wanita, dengan selisih sekitar 50%. Standar eror RD yang kecil menunjukkan bahwa selisih antara jenis kelamin yang terlihat memang benar ada di populasi, sehingga estimasi tersebut dapat dikatakan cukup akurat. Uji Z yang signifikan artinya terdapat beda peluang antara jenis kelamin laki-laki dan perempuan terhadap preferensi produk A.
Uji asosiasi resiko relatif (RR): Berdasarkan hasil perhitungan, nilai RR menunjukkan bahwapria memiliki kemungkinan 2,5 kali lebih besar untuk menyukai produk A dibandingkan wanita. Standar eror RR yang kecil menunjukkan bahwa resiko relatif antara jenis kelamin yang terlihat memang benar ada di populasi, sehingga estimasi tersebut dapat dikatakan cukup akurat. Uji Z yang signifikan artinya terdapat perbedaan yang signifikan dalam risiko menyukai produk A antara pria dan wanita.
Uji asosiasi odds ratio (OR): Berdasarkan hasil perhitungan, nilai OR menunjukkan bahwapeluang pria menyukai produk A adalah 10 kali lebih besar dibanding wanita. Standar eror OR yang kecil menunjukkan bahwa OR antara jenis kelamin yang terlihat memang benar ada di populasi, sehingga estimasi tersebut dapat dikatakan cukup akurat. Uji Z yang signifikan artinya terdapat asosiasi antara jenis kelamin terhadap preferensi produk A.
Partisi chi-square memungkinkan kita untuk memecah statistik chi-square ke dalam komponen-komponen yang lebih kecil, terutama untuk tabel yang lebih besar. Namun, untuk tabel 2x2, partisi ini tidak relevan karena derajat kebebasannya hanya 1, sehingga statistik chi-square tidak dapat dipartisi lebih lanjut.
Contoh soal:
Sebuah survei dilakukan untuk mengetahui apakah terdapat hubungan antara jenis kelamin (Pria/Wanita) dan preferensi grup band (girl grup korea, boy grup korea, dan JKT48).
Berikut data yang diperoleh dari 600 responden:
| Gender | Girl Grup Korea | Boy Grup Korea | JKT 48 | Total |
|---|---|---|---|---|
| Pria | 80 | 20 | 200 | 300 |
| Wanita | 90 | 130 | 80 | 300 |
| Total | 170 | 150 | 280 | 600 |
Tentukan:
1. Uji chi-square untuk mengetahui hubungan antara jenis kelamin dan preferensi grup band.
2. Lakukan uji partisi chi-square:
- Partisi 1: girl grup korea vs boy grup korea
- Partisi 2: grup band korea (girl grup korea + boy grup korea) vs grup band indonesia (JKT 48)
Uji Chi-Square (Induk)
Hipotesis
Tingkat Signifikansi
alfa = 5% = 0,05
Perhitungan R
# Membuat tabel kontingensi
tabel <- matrix(c(80, 20, 200,
90, 130, 80),
nrow = 2, byrow = TRUE)
colnames(tabel) <- c("girl grup korea", "boy grup korea", "JKT 48")
rownames(tabel) <- c("Pria", "Wanita")
tabel <- as.table(tabel)
# Menampilkan tabel
tabel
## girl grup korea boy grup korea JKT 48
## Pria 80 20 200
## Wanita 90 130 80
# Uji chi-square
uji_chi <- chisq.test(tabel)
uji_chi
##
## Pearson's Chi-squared test
##
## data: tabel
## X-squared = 132.68, df = 2, p-value < 2.2e-16
Kriteria Uji
Tolak H0 jika nilai p-value < alfa
Keputusan
Karena p-value (2,2e-16) < alfa (0,05), maka tolak H0. Artinya, terdapat hubungan antara variabel jenis kelamin terhadap preferensi grup band.
Uji Partisi (Setiap 2x2)
Hipotesis
Hipotesis untuk setiap partisi:
tabel1 <- matrix(c(80, 20,
90, 130),
nrow = 2, byrow = TRUE)
colnames(tabel1) <- c("girl grup korea", "boy grup korea")
rownames(tabel1) <- c("Pria", "Wanita")
tabel1 <- as.table(tabel1)
chisq.test(tabel1)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabel1
## X-squared = 40.632, df = 1, p-value = 1.837e-10
tabel2 <- matrix(c(100, 200,
220, 80),
nrow = 2, byrow = TRUE)
colnames(tabel2) <- c("girl grup korea + boy grup korea", "JKT 48")
rownames(tabel2) <- c("Pria", "Wanita")
tabel2 <- as.table(tabel2)
chisq.test(tabel2)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabel2
## X-squared = 94.828, df = 1, p-value < 2.2e-16
Kriteria Uji
Tolak H0 jika p-value < alfa
Keputusan
Karena p-value partisi pertama (1,837e-10) < alfa (0,05), maka tolak H0. Artinya, pria dan wanita punya preferensi berbeda antara girl grup korea dan boy grup korea.
Karena p-value partisi kedua (2,2e-16)< alfa (0,05), maka tolak H0. Artinya, pria dan wanita punya preferensi berbeda antara grup band korea dan grup band indonesia.
Kesimpulan
Uji Likelihood Ratio digunakan sebagai alternatif uji Chi-Square
untuk tabel besar.
\[ G^2 = 2 \sum O_{ij}
\ln\left(\frac{O_{ij}}{E_{ij}}\right) \]
Contoh soal:
Sebuah perusahaan melakukan survei terhadap 1000 orang untuk mengetahui apakah ada hubungan antara status merokok (Perokok / Bukan Perokok) dan kejadian penyakit jantung (Ya / Tidak).
Hasil survei tersebut disajikan dalam tabel berikut:
| Penyakit Jantung (Ya) | Tidak (Tidak) | Total | |
|---|---|---|---|
| Perokok | 180 | 320 | 500 |
| Bukan Perokok | 90 | 410 | 500 |
| Total | 270 | 730 | 1000 |
Apakah ada hubungan (tidak independen) antara status merokok dan kejadian penyakit jantung menggunakan uji likelihood ratio G².
Hipotesis
Perhitungan R
# Matriks data observasi
observed <- matrix(c(180, 320,
90, 410),
nrow = 2, byrow = TRUE)
colnames(observed) <- c("Penyakit_Jantung", "Tidak")
rownames(observed) <- c("Perokok", "Bukan_Perokok")
# Frekuensi total
total <- sum(observed)
# Frekuensi baris dan kolom
row_totals <- rowSums(observed)
col_totals <- colSums(observed)
# Matriks expected frequency
dataexpected <- chisq.test(observed)$expected
# G^2 manual
G2 <- 2 * sum(observed * log(observed / dataexpected))
critical_value <- qchisq(0.95, df = 1)
print(G2)
## [1] 41.706
print(critical_value)
## [1] 3.841459
Kriteria Uji
Tolak H0 jika nilai G² > nilai kritis χ²
Keputusan
Karena G² (41,706) > nilai kritis χ² (3,841), maka tolak H0.
Kesimpulan
Berdasarkan hasil Uji Likelihood Ratio di atas, dapat disimpulkan bahwa terdapat hubungan antara variabel perokok dan penyakit jantung.
Uji Fisher Exact digunakan untuk sampel kecil atau sel dengan
frekuensi <5.
\[ P =
\frac{(a+b)!(c+d)!(a+c)!(b+d)!}{n!a!b!c!d!} \]
Distribusi Geometrik
Distribusi hipergeometrik mendasari uji Exact Fisher. Distribusi ini digunakan untuk menghitung probabilitas dalam tabel kontingensi 2x2 ketika margin (total baris dan kolom) tetap. Probabilitas untuk mendapatkan \(n_{11}\) adalah:
\[ P(n_{11}) = \frac{\binom{n_{1.}}{n_{11}} \binom{n_{2.}}{n_{.1} - n_{11}}}{\binom{n}{n_{.1}}} \]
Sebuah penelitian kecil dilakukan untuk mengetahui apakah ada hubungan antara jenis kelamin dan keberhasilan lulus ujian matematika pada sekelompok siswa di kelas.
Berikut adalah hasil pengamatan terhadap total 12 siswa:
| Lulus | Tidak Lulus | Total | |
|---|---|---|---|
| Pria | 3 | 2 | 5 |
| Wanita | 1 | 6 | 7 |
| Total | 4 | 8 | 12 |
Apakah ada hubungan antara jenis kelamin dan kelulusan ujian matematika menggunakan uji Fisher’s Exact?
Hipotesis
Taraf Signifikansi alfa = 5% = 0,05
Perhitungan Manual (Distribusi Hipergeometrik)
Menggunakan distribusi hipergeometrik:
\[ P = \frac{\binom{a+b}{a} \binom{c+d}{c}}{\binom{n}{a+c}} \]
Dengan: - \(a = 3\), \(b = 2\), \(c = 1\), \(d = 6\) - \(n = 12\)
\[ P = \frac{\binom{5}{3} \cdot \binom{7}{1}}{\binom{12}{4}} = \frac{10 \cdot 7}{495} = \frac{70}{495} \approx 0.1414 \]
Perhitungan Manual (Distribusi Hipergeometrik)
Kita akan menghitung probabilitas untuk semua tabel 2x2 yang memiliki marginal total yang sama, yaitu: - Total pria: 5 - Total wanita: 7 - Total lulus: 4 - Total tidak lulus: 8
# Tabel empiris: Pria Lulus = 3, Wanita Lulus = 1
# Probabilitas = P(X = 3)
choose(5, 3) * choose(7, 1) / choose(12, 4) # Tabel empiris
## [1] 0.1414141
# Semua kombinasi yang mungkin: Pria lulus = 0, 1, 2, 3, 4 (max 4 karena total lulus 4)
probs <- c(
choose(5, 0) * choose(7, 4), # 0 pria lulus
choose(5, 1) * choose(7, 3), # 1
choose(5, 2) * choose(7, 2), # 2
choose(5, 3) * choose(7, 1), # 3 (empiris)
choose(5, 4) * choose(7, 0) # 4
) / choose(12, 4)
probs
## [1] 0.07070707 0.35353535 0.42424242 0.14141414 0.01010101
p_value_manual <- sum(probs[probs <= probs[4]]) # probs[4] adalah P(X = 3)
p_value_manual
## [1] 0.2222222
Perhitungan dengan R (Package fisher.test)
# Matriks data observasi
data <- matrix(c(3, 2,
1, 6),
nrow = 2, byrow = TRUE)
colnames(data) <- c("Lulus", "Tidak Lulus")
rownames(data) <- c("Pria", "Wanita")
data
## Lulus Tidak Lulus
## Pria 3 2
## Wanita 1 6
fisher.test(data)
##
## Fisher's Exact Test for Count Data
##
## data: data
## p-value = 0.2222
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.3527639 546.2288490
## sample estimates:
## odds ratio
## 7.218098
Kriteria uji
Tolak H0 jika p-value < alfa
Keputusan
Karena p-value (0,22) > alfa (0,05) maka terima H0
Kesimpulan
Berdasarkan hasil Uji Exact Fisher di atas, dapat disimpulkan bahwa tidak terdapat hubungan antara variabel jenis kelamin dan keberhasilan lulus ujian matematika.
Residual Pearson digunakan untuk mengukur penyimpangan observasi dari
ekspektasi.
\[ e_{ij} = \frac{O_{ij} -
E_{ij}}{\sqrt{E_{ij}}} \]
Residual terstandarisasi lebih memperhitungkan variabilitas
marginal.
\[ r_{ij} = \frac{O_{ij} -
E_{ij}}{\sqrt{E_{ij}(1 - p_{i+})(1 - p_{+j})}} \]
Deteksi Outlier:
- \(|r_{ij}| > 3\) → outlier
signifikan.
Contoh soal:
Sebuah survei dilakukan untuk mengetahui hubungan antara jenis kelamin dan preferensi terhadap minuman berenergi pada 200 responden.
| Suka | Tidak Suka | Jumlah | |
|---|---|---|---|
| Laki-laki | 95 | 15 | 110 |
| Perempuan | 35 | 55 | 90 |
| Jumlah | 130 | 70 | 200 |
Lakukan analisis residual Pearson dan Standardized (Adjusted) untuk mendeteksi kemungkinan outlier!
Hipotesis
Perhitungan R
# Matriks Observasi
observed <- matrix(c(30, 5, 40, 30), nrow = 2, byrow = TRUE)
colnames(observed) <- c("Suka", "Tidak Suka")
rownames(observed) <- c("Laki-laki", "Perempuan")
#Perhitungan manual
# Total baris dan kolom
baris_total <- rowSums(data)
kolom_total <- colSums(data)
total <- sum(data)
# Matriks ekspektasi (nilai harapan)
ekspektasi <- outer(baris_total, kolom_total) / total
ekspektasi
## Lulus Tidak Lulus
## Pria 1.666667 3.333333
## Wanita 2.333333 4.666667
# Pearson residual: (observed - expected) / sqrt(expected)
residual_pearson <- (data - ekspektasi) / sqrt(ekspektasi)
residual_pearson
## Lulus Tidak Lulus
## Pria 1.0327956 -0.7302967
## Wanita -0.8728716 0.6172134
#Perhitungan dengan package
# Uji Chi-Square
chisq.test(observed)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: observed
## X-squared = 7.3339, df = 1, p-value = 0.006767
# Hitung nilai harapan (expected)
expected <- chisq.test(observed)$expected
# Pearson Residuals
pearson_residual <- (observed - expected) / sqrt(expected)
pearson_residual
## Suka Tidak Suka
## Laki-laki 1.3801311 -1.951800
## Perempuan -0.9759001 1.380131
# Standardized Residuals (Adjusted)
n <- sum(observed)
row_prop <- rowSums(observed) /n
col_prop <- colSums(observed) /n
standardized_residual <- (observed - expected) /
sqrt(expected * (1 - row_prop) * (1 - col_prop))
standardized_residual
## Suka Tidak Suka
## Laki-laki 2.927700 -4.140393
## Perempuan -2.070197 2.927700
Kesimpulan
Berdasarkan perhitungan pearson residual, dapat disimpulkan bahwa sel laki-laki tidak suka memiliki frekuensi observasi yang jauh lebih rendah dari yang diharapkan. Sebaliknya, sel laki-laki suka dan sel perempuan tidak suka memiliki frekuensi observasi lebih tinggi dari yang diharapkan.
Berdasarkan perhitungan standardized residual, dapat disimpulkan bahwa terdapat sel outlier yaitu pada sel laki-laki tidak suka. Hal tersebut karena | r | > 3, sehingga sel tersebut dapat dianggap sebagai outlier signifikan.
Tabel kontingensi tiga arah digunakan untuk menganalisis hubungan antara tiga variabel kategorik, misalnya variabel X, Y, dan Z. Dalam analisis ini, kita dapat mempelajari hubungan antar variabel dengan mempertimbangkan efek variabel ketiga sebagai kontrol.
Misalkan:
Tabel kontingensi tiga arah dapat ditulis sebagai:
| \(Z = 1\) (muda) | \(Y=1\) (sakit) | \(Y=2\) (sehat) |
|---|---|---|
| \(X=1\) (terpapar) | \(n_{111}\) | \(n_{112}\) |
| \(X=2\) (tidak terpapar) | \(n_{121}\) | \(n_{122}\) |
| \(Z = 2\) (tua) | \(Y=1\) (sakit) | \(Y=2\) (sehat) |
|---|---|---|
| \(X=1\) (terpapar) | \(n_{211}\) | \(n_{212}\) |
| \(X=2\) (tidak terpapar) | \(n_{221}\) | \(n_{222}\) |
Di mana \(n_{ijk}\) adalah frekuensi pada \(X=i\), \(Y=j\), dan \(Z=k\).
Tabel parsial adalah tabel 2x2 yang dibuat dengan mengkondisikan pada level tertentu dari variabel ketiga (\(Z\)). Dalam contoh di atas, kita memiliki dua tabel parsial: satu untuk \(Z=1\) (muda) dan satu untuk \(Z=2\) (tua). Tabel parsial memungkinkan analisis hubungan antara \(X\) dan \(Y\) pada setiap level \(Z\), sehingga kita dapat melihat apakah hubungan tersebut dipengaruhi oleh \(Z\).
Tabel marjinal diperoleh dengan menjumlahkan frekuensi melintasi level variabel ketiga (\(Z\)), sehingga menghasilkan tabel 2x2 untuk \(X\) dan \(Y\) tanpa mempedulikan \(Z\). Frekuensi marjinal dihitung sebagai:
\[ n_{ij.} = {\sum_{k} n_{ijk}} \]
Contoh:
| \(Y=1\) (sakit) | \(Y=2\) (sehat) | |
|---|---|---|
| \(X=1\) (terpapar) | \(n_{11.} = n_{111} + n_{211}\) | \(n_{12.} = n_{112} + n_{212}\) |
| \(X=2\) (tidak terpapar) | \(n_{21.} = n_{121} + n_{221}\) | \(n_{22.} = n_{122} + n_{222}\) |
Tabel marjinal mengabaikan efek \(Z\) dan hanya menunjukkan hubungan langsung antara \(X\) dan \(Y\).
Distribusi peluang dalam tabel kontingensi tiga arah dapat diuraikan menjadi tiga jenis: joint probability, marginal probability, dan conditional probability.
Joint probability adalah peluang kejadian gabungan dari ketiga variabel \(X\), \(Y\), dan \(Z\):
\[ P(X=i, Y=j, Z=k) = \frac{n_{ijk}}{n} \]
Di mana \(n\) adalah total observasi.
Marginal probability adalah peluang suatu variabel tanpa mempedulikan variabel lain. Misalnya, peluang \(X=1\):
\[ P(X=1) = \frac{\sum_{j,k} n_{1jk}}{n} \]
Conditional probability adalah peluang suatu kejadian diberikan kejadian lain. Misalnya, peluang \(Y=1\) diberikan \(X=1\) dan \(Z=1\):
\[ P(Y=1 | X=1, Z=1) = \frac{n_{111}}{n_{11.}} \]
Di mana \(n_{11.} = n_{111} + n_{112}\).
Ukuran asosiasi seperti Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR) dapat dihitung untuk setiap tabel parsial.
RD mengukur selisih risiko antara dua kelompok pada level tertentu \(Z\):
\[ \hat{\text{RD}}k = {\hat{p}_{1|1,k}} - {\hat{p}_{1|2,k}} \]
Di mana:
\[ \hat{p}_{1|1,k} = \frac{n_{k11}}{n_{k1.}}, \quad \hat{p}_{1|2,k} = \frac{n_{k21}}{n_{k2.}} \]
RR mengukur rasio risiko antara dua kelompok pada level tertentu \(Z\):
\[ \hat{\text{RR}}_k = \frac{\hat{p}_{1|1,k}}{\hat{p}_{1|2,k}} \]
OR mengukur rasio odds antara dua kelompok pada level tertentu \(Z\):
\[ \hat{\text{OR}}k = \frac{n_{k11} n_{k22}}{n_{k12} n_{k21}} \]
Ukuran asosiasi ini dapat dibandingkan antar level \(Z\) untuk melihat apakah hubungan antara \(X\) dan \(Y\) dipengaruhi oleh \(Z\).
Conditional independence terjadi ketika hubungan antara \(X\) dan \(Y\) tidak dipengaruhi oleh \(Z\), artinya \(X\) dan \(Y\) independen pada setiap level \(Z\). Secara formal, \(X\) dan \(Y\) dikatakan conditionally independent diberikan \(Z\) jika:
\[ P(Y=j | X=i, Z=k) = P(Y=j | Z=k) \quad \text{untuk semua } i, j, k \]
Dalam hal ini, ukuran asosiasi seperti OR akan konstan (atau mendekati konstan) di semua level $Z$. Untuk menguji conditional independence, kita dapat membandingkan OR antar tabel parsial atau menggunakan uji statistik seperti uji chi-square pada setiap tabel parsial.
Contoh soal:
Sebuah studi melakukan penelitian tentang kebiasaan merokok (X), usia responden (Z), dan penyakit jantung (Y). Data berikut diperoleh dari 150 responden yang dikelompokkan berdasarkan ketiga variabel tersebut:
| Usia (Z) | Merokok (X) | Penyakit Jantung (Y) | Ya | Tidak | Jumlah |
|---|---|---|---|---|---|
| Muda | Ya | 5 | 25 | 30 | |
| Muda | Tidak | 2 | 28 | 30 | |
| Dewasa | Ya | 10 | 20 | 30 | |
| Dewasa | Tidak | 4 | 26 | 30 | |
| Tua | Ya | 15 | 5 | 20 | |
| Tua | Tidak | 6 | 4 | 10 |
Tentukan peluang bersama, peluang marjinal, peluang bersyarat, dan ukuran asosiasinya
Usia Muda (\(Z=1\)):
\(P(X=1, Y=1, Z=1)\) (Merokok, Penyakit Jantung, Muda):
\[ P(X=1, Y=1, Z=1) = \frac{5}{150} = 0.0333 \]
\(P(X=1, Y=0, Z=1)\) (Merokok, Tidak Penyakit Jantung, Muda):
\[ P(X=1, Y=0, Z=1) = \frac{25}{150} = 0.1667 \]
\(P(X=0, Y=1, Z=1)\) (Tidak Merokok, Penyakit Jantung, Muda):
\[ P(X=0, Y=1, Z=1) = \frac{2}{150} = 0.0133 \]
\(P(X=0, Y=0, Z=1)\) (Tidak Merokok, Tidak Penyakit Jantung, Muda):
\[ P(X=0, Y=0, Z=1) = \frac{28}{150} = 0.1867 \]
Usia Dewasa (\(Z=2\)):
\(P(X=1, Y=1, Z=2)\):
\[ P(X=1, Y=1, Z=2) = \frac{10}{150} = 0.0667 \]
\(P(X=1, Y=0, Z=2)\):
\[ P(X=1, Y=0, Z=2) = \frac{20}{150} = 0.1333 \]
\(P(X=0, Y=1, Z=2)\):
\[ P(X=0, Y=1, Z=2) = \frac{4}{150} = 0.0267 \]
\(P(X=0, Y=0, Z=2)\):
\[ P(X=0, Y=0, Z=2) = \frac{26}{150} = 0.1733 \]
Usia Tua (\(Z=3\)):
\(P(X=1, Y=1, Z=3)\):
\[ P(X=1, Y=1, Z=3) = \frac{15}{150} = 0.1 \]
\(P(X=1, Y=0, Z=3)\):
\[ P(X=1, Y=0, Z=3) = \frac{5}{150} = 0.0333 \]
\(P(X=0, Y=1, Z=3)\):
\[ P(X=0, Y=1, Z=3) = \frac{6}{150} = 0.04 \]
\(P(X=0, Y=0, Z=3)\):
\[ P(X=0, Y=0, Z=3) = \frac{4}{150} = 0.0267 \]
data <- array(c(5, 2, 25, 28, 10, 4, 20, 26, 15, 6, 5, 4),
dim = c(2, 2, 3),
dimnames = list(
Merokok = c("Ya", "Tidak"),
PJ = c("Ya", "Tidak"),
Usia = c("Muda", "Dewasa", "Tua")))
# Hitung probabilitas bersama
total <- sum(data)
joint_prob <- data / total
ftable(joint_prob)
## Usia Muda Dewasa Tua
## Merokok PJ
## Ya Ya 0.03333333 0.06666667 0.10000000
## Tidak 0.16666667 0.13333333 0.03333333
## Tidak Ya 0.01333333 0.02666667 0.04000000
## Tidak 0.18666667 0.17333333 0.02666667
\(P(Z=1)\) (Usia Muda):
\[ P(Z=1) = \frac{30 + 30}{150} = \frac{60}{150} = 0.4 \]
\(P(Z=2)\) (Usia Dewasa):
\[ P(Z=2) = \frac{30 + 30}{150} = \frac{60}{150} = 0.4 \]
\(P(Z=3)\) (Usia Tua):
\[ P(Z=3) = \frac{20 + 10}{150} = \frac{30}{150} = 0.2 \]
\(P(X=1)\) (Merokok):
\[ P(X=1) = \frac{30 + 30 + 20}{150} = \frac{80}{150} = 0.5333 \]
\(P(X=0)\) (Tidak Merokok):
\[ P(X=0) = \frac{30 + 30 + 10}{150} = \frac{70}{150} = 0.4667 \]
\(P(Y=1)\) (Penyakit Jantung):
\[ P(Y=1) = \frac{5 + 2 + 10 + 4 + 15 + 6}{150} = \frac{42}{150} = 0.28 \]
\(P(Y=0)\) (Tidak Penyakit Jantung):
\[ P(Y=0) = \frac{25 + 28 + 20 + 26 + 5 + 4}{150} = \frac{108}{150} = 0.72 \]
# Hitung probabilitas bersama
total <- sum(data)
joint_prob <- data / total
# Hitung probabilitas marginal
marginal_X <- apply(joint_prob, 1, sum)
marjinal_pj <- margin.table(data, 2) / total
marginal_Z <- apply(joint_prob, 3, sum)
# Tampilkan hasil
marginal_X
## Ya Tidak
## 0.5333333 0.4666667
marjinal_pj
## PJ
## Ya Tidak
## 0.28 0.72
marginal_Z
## Muda Dewasa Tua
## 0.4 0.4 0.2
Usia Muda (\(Z=1\))
\(P(Y=1 | X=1, Z=1)\) (Penyakit Jantung diberikan Merokok, Usia Muda):
\[ P(Y=1 | X=1, Z=1) = \frac{5}{30} = 0.1667 \]
\(P(Y=1 | X=0, Z=1)\) (Penyakit Jantung diberikan Tidak Merokok, Usia Muda):
\[ P(Y=1 | X=0, Z=1) = \frac{2}{30} = 0.0667 \]
Usia Dewasa (\(Z=2\))
\(P(Y=1 | X=1, Z=2)\):
\[ P(Y=1 | X=1, Z=2) = \frac{10}{30} = 0.3333 \]
\(P(Y=1 | X=0, Z=2)\):
\[ P(Y=1 | X=0, Z=2) = \frac{4}{30} = 0.1333 \]
Usia Tua (\(Z=3\))
\(P(Y=1 | X=1, Z=3)\):
\[ P(Y=1 | X=1, Z=3) = \frac{15}{20} = 0.75 \]
\(P(Y=1 | X=0, Z=3)\):
\[ P(Y=1 | X=0, Z=3) = \frac{6}{10} = 0.6 \]
# Inisialisasi array untuk simpan hasil peluang bersyarat
cond_probs <- array(0, dim = c(2, 2, 3),
dimnames = list(
Merokok = c("Ya", "Tidak"),
PJ = c("Ya", "Tidak"),
Usia = c("Muda", "Dewasa", "Tua")
))
# Loop untuk setiap usia
for (usia in dimnames(data)$Usia) {
# Ambil sub-tabel untuk usia tersebut (2x2)
sub_tabel <- data[ , , usia]
# Hitung total responden pada usia tersebut
total_usia <- sum(sub_tabel)
# Hitung peluang bersyarat: P(Merokok, PJ | Usia)
cond_probs[ , , usia] <- sub_tabel / total_usia
}
# Cetak hasil peluang bersyarat untuk masing-masing usia
cat("Peluang Bersyarat untuk Usia Muda:\n")
## Peluang Bersyarat untuk Usia Muda:
print(round(cond_probs[ , , "Muda"], 3))
## PJ
## Merokok Ya Tidak
## Ya 0.083 0.417
## Tidak 0.033 0.467
cat("\nPeluang Bersyarat untuk Usia Dewasa:\n")
##
## Peluang Bersyarat untuk Usia Dewasa:
print(round(cond_probs[ , , "Dewasa"], 3))
## PJ
## Merokok Ya Tidak
## Ya 0.167 0.333
## Tidak 0.067 0.433
cat("\nPeluang Bersyarat untuk Usia Tua:\n")
##
## Peluang Bersyarat untuk Usia Tua:
print(round(cond_probs[ , , "Tua"], 3))
## PJ
## Merokok Ya Tidak
## Ya 0.5 0.167
## Tidak 0.2 0.133
Kesimpulan
Ukuran Asosasi
Usia Muda (\(Z=1\))
Risk Difference (RD):
\[ \hat{\text{RD}}{1} = \hat{p}{1|1,1} - \hat{p}*{1|2,1} = 0.1667 - 0.0667 = 0.1 \]
Relative Risk (RR):
\[ \hat{\text{RR}}{1} = \frac{\hat{p}{1|1,1}}{\hat{p}*{1|2,1}} = \frac{0.1667}{0.0667} = 2.5 \]
Odds Ratio (OR):
\[ \hat{\text{OR}}{1} = \frac{n{111} n*{122}}{n*{112} n*{121}} = \frac{5 \times 28}{25 \times 2} = \frac{140}{50} = 2.8 \]
Usia Dewasa (\(Z=2\))
Risk Difference (RD):
\[ \hat{\text{RD}}{2} = \hat{p}{1|1,2} - \hat{p}*{1|2,2} = 0.3333 - 0.1333 = 0.2 \]
Relative Risk (RR):
\[ \hat{\text{RR}}{2} = \frac{\hat{p}{1|1,2}}{\hat{p}*{1|2,2}} = \frac{0.3333}{0.1333} = 2.5 \]
Odds Ratio (OR):
\[ \hat{\text{OR}}{2} = \frac{n{211} n*{222}}{n*{212} n*{221}} = \frac{10 \times 26}{20 \times 4} = \frac{260}{80} = 3.25 \]
Usia Tua (\(Z=3\))
Risk Difference (RD):
\[ \hat{\text{RD}}{3} = \hat{p}{1|1,3} - \hat{p}*{1|2,3} = 0.75 - 0.6 = 0.15 \]
Relative Risk (RR):
\[ \hat{\text{RR}}{3} = \frac{\hat{p}{1|1,3}}{\hat{p}*{1|2,3}} = \frac{0.75}{0.6} = 1.25 \]
Odds Ratio (OR):
\[ \hat{\text{OR}}{3} = \frac{n{311} n*{322}}{n*{312} n*{321}} = \frac{15 \times 4}{5 \times 6} = \frac{60}{30} = 2 \]
# Data array
data <- array(c(5, 2, 25, 28, 10, 4, 20, 26, 15, 6, 5, 4),
dim = c(2, 2, 3),
dimnames = list(
Merokok = c("Ya", "Tidak"),
PJ = c("Ya", "Tidak"),
Usia = c("Muda", "Dewasa", "Tua")
))
# Buat fungsi untuk menghitung RD, RR, OR manual
hitung_ukuran_asosiasi <- function(tab) {
# Susun dalam format:
# a = Merokok = Ya, PJ = Ya
# b = Merokok = Ya, PJ = Tidak
# c = Merokok = Tidak, PJ = Ya
# d = Merokok = Tidak, PJ = Tidak
a <- tab["Ya", "Ya"]
b <- tab["Ya", "Tidak"]
c <- tab["Tidak", "Ya"]
d <- tab["Tidak", "Tidak"]
# Risiko (probabilitas) PJ
risk_exposed <- a / (a + b) # Merokok
risk_unexposed <- c / (c + d) # Tidak merokok
# Odds
odds_exposed <- a / b
odds_unexposed <- c / d
# Hitung ukuran asosiasi
rd <- risk_exposed - risk_unexposed # Risk Difference
rr <- risk_exposed / risk_unexposed # Risk Ratio
or <- odds_exposed / odds_unexposed # Odds Ratio
# Output
return(list(RD = rd, RR = rr, OR = or))
}
# Hitung untuk tiap usia
usia_levels <- dimnames(data)$Usia
hasil_asosiasi <- list()
for (usia in usia_levels) {
tab <- data[ , , usia]
hasil_asosiasi[[usia]] <- hitung_ukuran_asosiasi(tab)
}
# Tampilkan hasil
hasil_asosiasi
## $Muda
## $Muda$RD
## [1] 0.1
##
## $Muda$RR
## [1] 2.5
##
## $Muda$OR
## [1] 2.8
##
##
## $Dewasa
## $Dewasa$RD
## [1] 0.2
##
## $Dewasa$RR
## [1] 2.5
##
## $Dewasa$OR
## [1] 3.25
##
##
## $Tua
## $Tua$RD
## [1] 0.15
##
## $Tua$RR
## [1] 1.25
##
## $Tua$OR
## [1] 2
Kesimpulan
Kesimpulan: perokok memiliki kemungkinan lebih besar untuk terkena penyakit jantung, terutama perokok di usia muda dan dewasa.
Contoh soal:
Sebuah studi ingin mengetahui hubungan antara kebiasaan merokok (X), usia responden (Z), dan penyakit jantung (Y). Datanya ditampilkan dalam syntax R berikut:
| Usia (\(Z\)) | Merokok (\(X\)) | Penyakit Jantung (\(Y\)) | Ya | Tidak | Jumlah |
|---|---|---|---|---|---|
| Muda | Ya | 35 | 20 | 55 | |
| Muda | Tidak | 15 | 40 | 55 | |
| Tua | Ya | 45 | 25 | 70 | |
| Tua | Tidak | 10 | 30 | 40 |
Cara 1: membandingkan Odds Ratio
Usia Muda (\(Z=1\))
\(n_{111}\) = 35 (Merokok Ya, Penyakit Jantung Ya)
\(n_{112}\) = 20 (Merokok Ya, Penyakit Jantung Tidak)
\(n_{121}\) = 15 (Merokok Tidak, Penyakit Jantung Ya)
\(n_{122}\) = 40 (Merokok Tidak, Penyakit Jantung Tidak)
\[ \hat{\text{OR}}*{1} = \frac{35 \times 40}{20 \times 15} = \frac{1400}{300} = 4.6667 \]
Usia Tua (\(Z=2\))
\(n_{211}\) = 45 (Merokok Ya, Penyakit Jantung Ya)
\(n_{212}\) = 25 (Merokok Ya, Penyakit Jantung Tidak)
\(n_{221}\) = 10 (Merokok Tidak, Penyakit Jantung Ya)
\(n_{222}\) = 30 (Merokok Tidak, Penyakit Jantung Tidak)
\[ \hat{\text{OR}}*{2} = \frac{45 \times 30}{25 \times 10} = \frac{1350}{250} = 5.4 \]
data <- array(c(
35, 15, 20, 40, # Usia Muda
45, 10, 25, 30 # Usia Tua
), dim = c(2, 2, 2),
dimnames = list(
Merokok = c("Ya", "Tidak"),
PJ = c("Ya", "Tidak"),
Usia = c("Muda", "Tua")
))
calculate_or <- function(table) {
a <- table[1, 1] # Non Perguruan Tinggi & Non Sejahtera
b <- table[1, 2] # Non Perguruan Tinggi & Sejahtera
c <- table[2, 1] # Perguruan Tinggi & Non Sejahtera
d <- table[2, 2] # Perguruan Tinggi & Sejahtera
or_value <- (a * d) / (b * c)
return(or_value)
}
# Tabel parsial untuk Rural
muda_table <- data[, , "Muda"]
or_muda <- calculate_or(muda_table)
# Tabel parsial untuk Urban
tua_table <- data[, , "Tua"]
or_tua <- calculate_or(tua_table)
# Tabel marjinal (menjumlahkan lokasi)
marginal_table <- apply(data, c(1, 2), sum)
or_marginal <- calculate_or(marginal_table)
# Menampilkan hasil
cat("Odds Ratio Muda:", or_muda, "\n")
## Odds Ratio Muda: 4.666667
cat("Odds Ratio Tua:", or_tua, "\n")
## Odds Ratio Tua: 5.4
cat("Odds Ratio Marginal:", or_marginal, "\n")
## Odds Ratio Marginal: 4.977778
Karena odds ratio untuk setiap tabel parsial tidak sama dengan 1, maka tidak memenuhi syarat independensi kondisional. Artinya, variabel merokok dan sakit jantung tidak independen secara kondisional.
Cara 2: membandingkan Peluang Bersyarat
# Hitung peluang bersama
total <- sum(data)
joint_prob <- prop.table(data)
# Buat fungsi untuk cek conditional independence
check_conditional_independence <- function(data, usia) {
# Ambil slice untuk satu usia
slice <- data[ , , usia]
total_z <- sum(slice)
# Hitung P(X|Z) dan P(Y|Z)
px_given_z <- margin.table(slice, 1) / total_z # Merokok | Usia
py_given_z <- margin.table(slice, 2) / total_z # PJ | Usia
# Hitung P(X,Y|Z)
pxy_given_z <- slice / total_z
# Hitung P(X|Z) * P(Y|Z)
expected <- outer(px_given_z, py_given_z)
cat("\n=== Usia:", usia, "===\n")
print(round(pxy_given_z, 3))
cat("Expected (P(X|Z)*P(Y|Z)):\n")
print(round(expected, 3))
cat("Selisih mutlak:\n")
print(round(abs(pxy_given_z - expected), 3))
}
# Cek untuk setiap usia
check_conditional_independence(data, "Muda")
##
## === Usia: Muda ===
## PJ
## Merokok Ya Tidak
## Ya 0.318 0.182
## Tidak 0.136 0.364
## Expected (P(X|Z)*P(Y|Z)):
## PJ
## Merokok Ya Tidak
## Ya 0.227 0.273
## Tidak 0.227 0.273
## Selisih mutlak:
## PJ
## Merokok Ya Tidak
## Ya 0.091 0.091
## Tidak 0.091 0.091
check_conditional_independence(data, "Tua")
##
## === Usia: Tua ===
## PJ
## Merokok Ya Tidak
## Ya 0.409 0.227
## Tidak 0.091 0.273
## Expected (P(X|Z)*P(Y|Z)):
## PJ
## Merokok Ya Tidak
## Ya 0.318 0.318
## Tidak 0.182 0.182
## Selisih mutlak:
## PJ
## Merokok Ya Tidak
## Ya 0.091 0.091
## Tidak 0.091 0.091
Karena 𝑃(𝑋,𝑌|𝑍) ≠ 𝑃(𝑋|𝑍)𝑃(𝑌|𝑍), maka, variabel merokok dan sakit jantung tidak independen secara kondisional.
Cara 3: uji Chi-Square
Hipotesis
Taraf Signifikansi
alfa = 5% = 0,05
Perhitungan manual
Usia Muda (\(Z=1\))
Hitung frekuensi yang diharapkan:
\(e_{11}\) (Merokok Ya, Penyakit Jantung Ya):
\[ e*{11} = \frac{n*{1.} \cdot n*{.1}}{n} = \frac{55 \times 50}{110} = \frac{2750}{110} = 25 \]
\(e_{12}\) (Merokok Ya, Penyakit Jantung Tidak):
\[ e*{12} = \frac{n*{1.} \cdot n*{.2}}{n} = \frac{55 \times 60}{110} = \frac{3300}{110} = 30 \]
\(e_{21}\) (Merokok Tidak, Penyakit Jantung Ya):
\[ e*{21} = \frac{n*{2.} \cdot n*{.1}}{n} = \frac{55 \times 50}{110} = \frac{2750}{110} = 25 \]
\(e_{22}\) (Merokok Tidak, Penyakit Jantung Tidak):
\[ e*{22} = \frac{n*{2.} \cdot n*{.2}}{n} = \frac{55 \times 60}{110} = \frac{3300}{110} = 30 \]
Hitung statistik Chi-Square:
\[ \chi^2 = \frac{(n*{11} - e*{11})^2}{e*{11}} + \frac{(n*{12} - e*{12})^2}{e*{12}} + \frac{(n*{21} - e*{21})^2}{e*{21}} + \frac{(n*{22} - e*{22})^2}{e*{22}} \]
\[ \chi^2 = \frac{(35 - 25)^2}{25} + \frac{(20 - 30)^2}{30} + \frac{(15 - 25)^2}{25} + \frac{(40 - 30)^2}{30} \]
\[ \chi^2 = \frac{(10)^2}{25} + \frac{(-10)^2}{30} + \frac{(-10)^2}{25} + \frac{(10)^2}{30} \]
\[ \chi^2 = \frac{100}{25} + \frac{100}{30} + \frac{100}{25} + \frac{100}{30} \]
\[ \chi^2 = 4 + 3.3333 + 4 + 3.3333 = 14.6666 \]
Usia Tua (\(Z=2\))
Hitung frekuensi yang diharapkan:
\(e_{11}\) (Merokok Ya, Penyakit Jantung Ya):
\[ e*{11} = \frac{n*{1.} \cdot n*{.1}}{n} = \frac{70 \times 55}{110} = \frac{3850}{110} = 35 \]
\(e_{12}\) (Merokok Ya, Penyakit Jantung Tidak):
\[ e*{12} = \frac{n*{1.} \cdot n*{.2}}{n} = \frac{70 \times 55}{110} = \frac{3850}{110} = 35 \]
\(e_{21}\) (Merokok Tidak, Penyakit Jantung Ya):
\[ e*{21} = \frac{n*{2.} \cdot n*{.1}}{n} = \frac{40 \times 55}{110} = \frac{2200}{110} = 20 \]
\(e_{22}\) (Merokok Tidak, Penyakit Jantung Tidak):
\[ e*{22} = \frac{n*{2.} \cdot n*{.2}}{n} = \frac{40 \times 55}{110} = \frac{2200}{110} = 20 \]
Hitung statistik Chi-Square:
\[ \chi^2 = \frac{(45 - 35)^2}{35} + \frac{(25 - 35)^2}{35} + \frac{(10 - 20)^2}{20} + \frac{(30 - 20)^2}{20} \]
\[ \chi^2 = \frac{(10)^2}{35} + \frac{(-10)^2}{35} + \frac{(-10)^2}{20} + \frac{(10)^2}{20} \]
\[ \chi^2 = \frac{100}{35} + \frac{100}{35} + \frac{100}{20} + \frac{100}{20} \]
\[ \chi^2 = 2.8571 + 2.8571 + 5 + 5 = 15.7142 \]
Perhitungan R
# Uji Chi-Square untuk masing-masing strata usia
chisq_muda <- chisq.test(data[,,"Muda"])
chisq_muda
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: data[, , "Muda"]
## X-squared = 13.237, df = 1, p-value = 0.0002745
chisq_tua <- chisq.test(data[,,"Tua"])
chisq_tua
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: data[, , "Tua"]
## X-squared = 14.182, df = 1, p-value = 0.0001659
Kriteria Uji
Tolak H0 jika p-value < alfa
Keputusan
Karena p-value untuk kelompok usia muda dan tua < alfa, maka tolak H0
Kesimpulan
Berdasarkan perhitungan di atas, dapat disimpulkan bahwa variabel merokok dan sakit jantung tidak independen secara kondisional. Artinya, merokok dan penyakit jantung saling berhubungan di tiap strata usia (muda dan tua)
Inferensi pada tabel kontingensi tiga arah melibatkan analisis hubungan antara dua variabel kategorik (misalnya \(X\) dan \(Y\)) dengan mengontrol efek variabel ketiga (\(Z\)). Dalam konteks ini, kita dapat menguji hubungan antar variabel menggunakan pendekatan statistik seperti uji Cochran-Mantel-Haenszel (CMH), menghitung odds ratio bersama, dan menguji homogenitas odds ratio antar strata dengan uji Breslow-Day. Tabel kontingensi tiga arah biasanya terdiri dari beberapa tabel 2x2 yang dikondisikan pada level \(Z\).
Misalkan:
Setiap level \(Z\) membentuk tabel parsial 2x2:
| \(Z = k\) | \(Y=1\) (sakit) | \(Y=2\) (sehat) |
|---|---|---|
| \(X=1\) (terpapar) | \(n\_{k11}\) | \(n\_{k12}\) |
| \(X=2\) (tidak terpapar) | \(n\_{k21}\) | \(n\_{k22}\) |
Uji Cochran-Mantel-Haenszel (CMH) digunakan untuk menguji hubungan antara \(X\) dan \(Y\) dengan mengontrol efek \(Z\), dengan asumsi odds ratio homogen antar strata \(Z\).
Statistik CMH dihitung sebagai:
\[ \text{CMH} = \frac{\left( \sum_{k} \left( n_{k11} - e_{k11} \right) \right)^2}{\sum_{k} v_{k11}} \]
Di mana:
\(e_{k11}\) adalah frekuensi yang diharapkan di sel \((1,1)\) pada strata \(k\) di bawah hipotesis nol:
\[ e_{k11} = \frac{(n_{k1.} \cdot n_{k.1})}{n_{k..}} \]
\(v_{k11}\) adalah varians dari \(n_{k11}\) di bawah hipotesis nol:
\[ v_{k11} = \frac{n_{k1.} \cdot n_{k2.} \cdot n_{k.1} \cdot n_{k.2}}{(n_{k..})^2 (n_{k..} - 1)} \]
\(n_{k..}\) adalah total observasi pada strata \(k\).
Statistik CMH mengikuti distribusi \(\chi^2\) dengan derajat kebebasan 1. Jika nilai CMH besar (p-value kecil), kita tolak \(H_0\), menunjukkan adanya asosiasi antara \(X\) dan \(Y\) setelah mengontrol \(Z\).
Odds ratio bersama (common odds ratio) menggabungkan odds ratio dari setiap strata \(Z\) ke dalam satu estimasi tunggal, dengan asumsi homogenitas odds ratio antar strata. Odds ratio untuk strata \(k\) adalah:
\[ \hat{\text{OR}}_k = \frac{n_{k11} n_{k22}}{n_{k12} n_{k21}} \]
Odds ratio bersama Cochran-Mantel-Haenszel dihitung sebagai:
\[ \hat{\text{OR}}{\text{CMH}} = \frac{\sum{k} \left( \frac{n_{k11} n_{k22}}{n_{k..}} \right)}{\sum_{k} \left( \frac{n_{k12} n_{k21}}{n_{k..}} \right)} \]
Odds ratio ini memberikan estimasi asosiasi antara \(X\) dan \(Y\) yang disesuaikan dengan efek \(Z\). Jika \(\hat{\text{OR}}_{\text{CMH}} \neq 1\), ada indikasi asosiasi antara \(X\) dan \(Y\) setelah mengontrol \(Z\).
Uji Breslow-Day digunakan untuk menguji homogenitas odds ratio antar strata \(Z\), yaitu apakah odds ratio \(\text{OR}_k\) sama untuk semua level \(Z\).
Hipotesis
\(H_0\): Odds ratio homogen antar strata (\(\text{OR}_1 = \text{OR}_2 = \cdots = \text{OR}_K\)).
\(H_1\): Odds ratio tidak homogen antar strata.
Statistik Breslow-Day dihitung sebagai:
\[ \text{BD} = \sum_{k} \frac{(n_{k11} - \hat{e}_{k11})^2}{\hat{v}_{k11}} \]
Di mana:
\(\hat{e}{k11}\) adalah frekuensi yang diharapkan di sel \((1,1)\) pada strata \(k\), dihitung menggunakan odds ratio bersama \(\hat{\text{OR}}{\text{CMH}}\) dan memecahkan persamaan berikut untuk \(n_{k11}\):
\[ \hat{\text{OR}}{\text{CMH}} = \frac{\hat{e}_{k11} (n_{k..} - n_{k1.} - n_{k.1} + \hat{e}_{k11})}{(n_{k1.} - \hat{e}_{k11})(n_{k.1} - \hat{e}_{k11})} \]
\(\hat{v}{k11}\) adalah varians dari \(\hat{e}{k11}\):
\[ \hat{v}{k11} = \left( \frac{1}{\hat{e}{k11}} + \frac{1}{n_{k1.} - \hat{e}{k11}} + \frac{1}{n{k.1} - \hat{e}{k11}} + \frac{1}{n{k..} - n_{k1.} - n_{k.1} + \hat{e}_{k11}} \right)^{-1} \]
Statistik BD mengikuti distribusi \(\chi^2\) dengan derajat kebebasan \(K-1\), di mana \(K\) adalah jumlah strata. Jika p-value kecil, kita tolak \(H_0\), menunjukkan bahwa odds ratio tidak homogen antar strata.
Contoh soal:
(Contoh soal sebelumnya) Sebuah studi ingin mengetahui hubungan antara kebiasaan merokok (X), usia responden (Z), dan penyakit jantung (Y). Datanya ditampilkan dalam syntax R berikut:
| Usia (\(Z\)) | Merokok (\(X\)) | Penyakit Jantung (\(Y\)) | Ya | Tidak | Jumlah |
|---|---|---|---|---|---|
| Muda | Ya | 35 | 20 | 55 | |
| Muda | Tidak | 15 | 40 | 55 | |
| Tua | Ya | 45 | 25 | 70 | |
| Tua | Tidak | 10 | 30 | 40 |
Hipotesis
H₀: Odds ratio antara merokok dan sakit jantung sama di setiap kategori usia (muda & tua)
H₁: Terdapat perbedaan odds ratio antara merokok dan sakit jantung di setiap kategori usia (muda & tua)
Taraf Signifikansi
alfa = 5% = 0,05
Statistik Uji
library(epitools)
library(DescTools)
## Warning: package 'DescTools' was built under R version 4.3.3
breslow_test <- BreslowDayTest(data)
print(breslow_test)
##
## Breslow-Day test on Homogeneity of Odds Ratios
##
## data: data
## X-squared = 0.058258, df = 1, p-value = 0.8093
breslowday.test<-function(x){
#Findthe commonORbasedonMantel-Haenszel
or.hat.mh <-mantelhaen.test(x)$estimate
#Numberofstrata
K <-dim(x)[3]
#Valueof theStatistic
X2.HBD<-0
#Valueof aj,tildeajandVar.aj
a <-tildea<-Var.a<-numeric(K)
for(j in 1:K){
#Findmarginalsoftablej
mj<-apply(x[,,j],MARGIN=1,sum)
nj<-apply(x[,,j],MARGIN=2,sum)
#Solve fortilde(a)_j
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))
#Taketheroot,whichfulfills0<tilde(a)_j<=min(n1_j,m1_j)
tildeaj<-sols[(0 < sols)& (sols<= min(nj[1],mj[1]))]
#Observedvalue
aj<-x[1,1,j]
#Determineotherexpectedcellentries
tildebj<-mj[1]-tildeaj
tildecj<-nj[1]-tildeaj
tildedj<-mj[2]-tildecj
#Compute\hat{\Var}(a_j|\widehat{\OR}_MH)
Var.aj<-(1/tildeaj+1/tildebj+ 1/tildecj+ 1/tildedj)^(-1)
#Computecontribution
X2.HBD<-X2.HBD + as.numeric((aj-tildeaj)^2 /Var.aj)
#Assignfoundvalueforlatercomputations
a[j]<-aj; tildea[j]<-tildeaj;Var.a[j]<-Var.aj
}
#ComputeTaronecorrectedtest
X2.HBDT<-as.numeric(X2.HBD-(sum(a)-sum(tildea))^2/sum(Var.aj))
#Computep-valuebasedontheTaronecorrectedtest
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)
}
print.bdtest<-function(x){
cat("BreslowandDaytest(withTaronecorrection):\n")
cat("Breslow-DayX-squared =",x$X2.HBD,"\n")
cat("Breslow-Day-TaroneX-squared =",x$X2.HBDT,"\n\n")
cat("TestfortestofacommonOR:p-value=",x$p,"\n\n")
}
breslowday.test(data)
## BreslowandDaytest(withTaronecorrection):
## Breslow-DayX-squared = 0.05825772
## Breslow-Day-TaroneX-squared = 0.05825269
##
## TestfortestofacommonOR:p-value= 0.8092793
Kriteria Uji
Tolak H0 jika p-value < alfa (0.05)
Keputusan
Karena p-value (0,8093) > alfa (0,05), maka terima H0
Kesimpulan
Berdasarkan perhitungan di atas, dapat disimpulkan bahwa tidak terdapat cukup bukti untuk menyatakan bahwa odds ratio untuk setiap strata berbeda. Artinya, efek merokok terhadap penyakit jantung dianggap sama, baik untuk usia muda maupun usia tua. Karena data tersebut memenuhi asumsi homogenitas, maka analisis dapat dilanjutkan ke Uji CMH.
Hipotesis
H₀: Tidak ada hubungan antara merokok dan penyakit jantung di strata usia.
H₁: Terdapat hubungan antara merokok dan penyakit jantung di strata usia.
Taraf Signifikansi
alfa = 5% = 0,05
Statistik Uji
Perhitungan manual
Untuk Usia Muda (\(k=1\)):
\[ e*{111} = \frac{(n*{111} + n*{112})(n*{111} + n*{121})}{n*{1}} = \frac{55 \times 50}{110} = 25 \]
\[ \text{Var}(n*{111}) = \frac{(n*{111} + n*{112})(n*{111} + n*{121})(n*{112} + n*{122})(n*{121} + n*{122})}{n*{1}^2 (n*{1} - 1)} \]
\[ \text{Var}(n*{111}) = \frac{55 \times 50 \times 60 \times 55}{110^2 \times 109} = \frac{9075000}{1321100} \approx 6.869 \]
\[ n*{111} - e*{111} = 35 - 25 = 10 \]
Untuk Usia Tua (\(k=2\)):
\[ e*{211} = \frac{(n*{211} + n*{212})(n*{211} + n*{221})}{n*{2}} = \frac{70 \times 55}{110} = 35 \]
\[ \text{Var}(n*{211}) = \frac{70 \times 55 \times 55 \times 40}{110^2 \times 109} = \frac{8470000}{1321100} \approx 6.411 \]
\[ n*{211} - e*{211} = 45 - 35 = 10 \]
Jumlahkan:
\[ \sum_{k} (n*{k11} - e*{k11}) = 10 + 10 = 20 \]
\[ \sum_{k} \text{Var}(n*{k11}) = 6.869 + 6.411 = 13.28 \]
\[ \chi^2_{\text{CMH}} = \frac{(20)^2}{13.28} = \frac{400}{13.28} \approx 30.12 \]
Perhitungan R
library(vcdExtra)
## Warning: package 'vcdExtra' was built under R version 4.3.3
## Loading required package: vcd
## Warning: package 'vcd' was built under R version 4.3.3
## Loading required package: grid
##
## Attaching package: 'vcd'
## The following object is masked from 'package:epitools':
##
## oddsratio
## Loading required package: gnm
## Warning: package 'gnm' was built under R version 4.3.3
##
## Attaching package: 'vcdExtra'
## The following object is masked from 'package:epitools':
##
## expand.table
## The following object is masked from 'package:dplyr':
##
## summarise
cmh_test <- mantelhaen.test(data)
print(cmh_test)
##
## Mantel-Haenszel chi-squared test with continuity correction
##
## data: data
## Mantel-Haenszel X-squared = 28.584, df = 1, p-value = 8.971e-08
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
## 2.769001 9.028528
## sample estimates:
## common odds ratio
## 5
# Inisialisasi nilai
numerator_sum <- 0
variance_sum <- 0
# Loop untuk setiap strata
for (k in 1:dim(data)[3]) {
a <- data["Ya", "Ya", k]
b <- data["Ya", "Tidak", k]
c <- data["Tidak", "Ya", k]
d <- data["Tidak", "Tidak", k]
n <- a + b + c + d
r1 <- a + b
r2 <- c + d
c1 <- a + c
c2 <- b + d
E_a <- r1 * c1 / n
Var_a <- (r1 * r2 * c1 * c2) / (n^2 * (n - 1))
numerator_sum <- numerator_sum + (a - E_a)
variance_sum <- variance_sum + Var_a
}
# Continuity correction
numerator_corrected <- numerator_sum - 0.5 * sign(numerator_sum)
# Statistik CMH
X2_cmh <- (numerator_corrected^2) / variance_sum
p_value <- pchisq(X2_cmh, df = 1, lower.tail = FALSE)
# Output
cat("CMH X-squared:", X2_cmh, "\n")
## CMH X-squared: 28.58431
cat("p-value:", p_value, "\n")
## p-value: 8.970611e-08
Kriteria Uji
Tolak H0 jika p-value < alfa (0.05)
Keputusan
Karena p-value (8,971e-08) < alfa(0,05), maka tolak H0
Kesimpulan
Berdasarkan perhitungan di atas, dapat disimpulkan bahwa terdapat hubungan yang signifikan antara merokok dan penyakit jantung setelah memperhitungkan strata usia (muda dan tua).
Perhitungan manual
\[ \hat{\text{OR}}_{\text{CMH}} = \frac{\sum_{k} \frac{n*{k11} n*{k22}}{n*{k}}}{\sum_{k} \frac{n*{k12} n*{k21}}{n*{k}}} = 5 \]
# Inisialisasi
K <- dim(data)[3] # jumlah strata (2)
numerator <- 0
denominator <- 0
n_total <- rep(0, K)
# Variabel pendukung untuk SE log odds ratio
se1 <- 0
se2 <- 0
se3 <- 0
for (k in 1:K) {
# Ambil sel dari tabel 2x2 ke-k
n11 <- data[1,1,k]
n12 <- data[1,2,k]
n21 <- data[2,1,k]
n22 <- data[2,2,k]
n <- n11 + n12 + n21 + n22
n_total[k] <- n
# Numerator dan Denominator dari Odds Ratio MH
numerator <- numerator + (n11 * n22) / n
denominator <- denominator + (n12 * n21) / n
# Bagian SE
se1 <- se1 + ((n11 + n12) * (n11 * n22) * (n^2)) / (2 * (n11 + n12) * (n * (n12)))
se2 <- se2 + ((n11 + n22) * (n11 + n12) + (n12 + n21) * (n11 + n22)) * (n11 * n22) / (2 * (n * n) * (n12 * n21))
se3 <- se3 + ((n12 + n21) * (n12 * n21) * (n^2)) / (2 * (n12 * n21)^2)
}
# Odds Ratio MH
theta_MH <- numerator / denominator
# Log Odds Ratio dan SE-nya (manual rumus sesuai gambar)
log_theta <- log(theta_MH)
# Karena perhitungan SE di atas terlalu panjang secara rumus, kita bisa pakai pendekatan var(log(OR_MH)) manual:
var_log_theta <- 1 / numerator + 1 / denominator
se_log_theta <- sqrt(var_log_theta)
# Interval Kepercayaan 95%
z <- qnorm(0.975)
lower <- log_theta - z * se_log_theta
upper <- log_theta + z * se_log_theta
# Cetak hasil
cat("Odds Ratio MH:", round(theta_MH, 4), "\n")
## Odds Ratio MH: 5
cat("Standard Error log(OR_MH):", round(se_log_theta, 4), "\n")
## Standard Error log(OR_MH): 0.4899
cat("95% CI for log(OR_MH): [", round(lower, 4), ",", round(upper, 4), "]\n")
## 95% CI for log(OR_MH): [ 0.6493 , 2.5696 ]
cat("95% CI for OR_MH: [", round(exp(lower), 4), ",", round(exp(upper), 4), "]\n")
## 95% CI for OR_MH: [ 1.9141 , 13.0609 ]
Berdasarkan perhitungan di atas, dapat disimpulkan bahwa setelah dikontrol berdasarkan usia, seorang perokok memiliki risiko 5 kali lebih tinggi untuk terkena penyakit jantung dibandingkan non perokok. Standard error sebesar 0.4899 menunjukkan estimasi odds ratio yang cukup baik. Interval kepercayaan menunjukkan bahwa dengan keyakinan 95%, odds ratio bersama berada antara 1.91 hingga 13.06. Karena seluruh interval berada di atas 1, maka hubungan antara status penyakit jantung dan merokok adalah signifikan secara statistik.
Generalized Linear Model (GLM) adalah kerangka statistik yang memperluas model regresi linier untuk menangani variabel respons yang tidak mengikuti distribusi normal, seperti data kategorik, data hitung, atau data biner. GLM terdiri dari tiga komponen utama: distribusi probabilitas dari keluarga eksponensial, fungsi link, dan model linier untuk prediktor. GLM sangat berguna dalam analisis data kategorik, dengan dua model utama yang sering digunakan: regresi logistik untuk data biner dan regresi Poisson untuk data hitung.
Distribusi dari keluarga eksponensial (exponential family) adalah dasar dari GLM. Distribusi ini mencakup banyak distribusi umum seperti normal, binomial, dan Poisson. Bentuk umum dari distribusi eksponensial adalah:
\[ f(y; \theta, \phi) = \exp\left( \frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi) \right) \]
Di mana:
\(y\): Variabel respons.
\(\theta\): Parameter kanonik (natural parameter).
\(\phi\): Parameter dispersi (dispersion parameter).
\(b(\theta)\): Fungsi yang menentukan distribusi.
\(a(\phi)\): Fungsi skala.
\(c(y, \phi)\): Fungsi normalisasi.
Distribusi yang termasuk dalam keluarga eksponensial memiliki sifat-sifat statistik yang memudahkan estimasi parameter, seperti hubungan langsung antara ekspektasi dan varians:
Ekspektasi: \(E(Y) = \mu = b'(\theta)\)
Varians: \(\text{Var}(Y) = b''(\theta) \cdot a(\phi)\)
Contoh distribusi dalam keluarga eksponensial yang relevan untuk data kategorik adalah:
Binomial (untuk regresi logistik).
Poisson (untuk regresi Poisson).
Regresi logistik digunakan untuk memodelkan variabel respons biner (misalnya, 0/1 atau sukses/gagal). Dalam GLM, regresi logistik mengasumsikan bahwa variabel respons \(Y\) mengikuti distribusi binomial.
Model regresi logistik:
Misalkan \(Y_i \sim \text{Binomial}(1, \pi_i)\), di mana \(\pi_i\) adalah probabilitas sukses untuk observasi ke-\(i\). Model regresi logistik menghubungkan \(\pi_i\) dengan prediktor linier \(\mathbf{x}_i^\top \beta\) melalui fungsi link logit:
\[ \log\left( \frac{\pi_i}{1 - \pi_i} \right) = \mathbf{x}_i^\top \beta \]
Di mana:
\(\frac{\pi_i}{1 - \pi_i}\): Odds sukses.
\(\mathbf{x}_i^\top \beta\): Kombinasi linier dari prediktor \(\mathbf{x}_i\) dan parameter \(\beta\).
Dari persamaan di atas, probabilitas \(\pi_i\) dapat ditulis sebagai:
\[ \pi_i = \frac{\exp(\mathbf{x}_i^\top \beta)}{1 + \exp(\mathbf{x}_i^\top \beta)} \]
Fungsi Link
Fungsi link logit adalah:
\[ g(\pi_i) = \log\left( \frac{\pi_i}{1 - \pi_i} \right) \]
Fungsi ini memetakan probabilitas \(\pi_i \in (0,1)\) ke ruang bilangan real \((-\infty, \infty)\).
Regresi Poisson digunakan untuk memodelkan variabel respons berupa data hitung (count data), seperti jumlah kejadian dalam periode waktu tertentu. Dalam GLM, regresi Poisson mengasumsikan bahwa variabel respons \(Y\) mengikuti distribusi Poisson.
Model regresi poisson:
Misalkan \(Y_i \sim \text{Poisson}(\lambda_i)\), di mana \(\lambda_i\) adalah rata-rata kejadian untuk observasi ke-\(i\). Model regresi Poisson menghubungkan \(\lambda_i\) dengan prediktor linier melalui fungsi link log:
\[ \log(\lambda_i) = \mathbf{x}_i^\top \beta \]
Sehingga:
\[ \lambda_i = \exp(\mathbf{x}_i^\top \beta) \]
Fungsi Link
Fungsi link log adalah:
\[ g(\lambda_i) = \log(\lambda_i) \]
Fungsi ini memetakan rata-rata \(\lambda_i \in (0, \infty)\) ke ruang bilangan real \((-\infty, \infty)\).
Sifat Distribusi Poisson
Dalam distribusi Poisson, ekspektasi dan varians sama:
\[ E(Y_i) = \text{Var}(Y_i) = \lambda_i \]
Selain distribusi dari keluarga eksponensial dan fungsi link, GLM juga memiliki komponen ketiga, yaitu prediktor linier \(\mathbf{x}_i^\top \beta\), yang memungkinkan fleksibilitas dalam memodelkan hubungan antara variabel respons dan prediktor. GLM memungkinkan analisis data kategorik dengan cara yang lebih fleksibel dibandingkan regresi linier biasa, terutama karena:
Tidak mensyaratkan asumsi normalitas pada variabel respons.
Mengakomodasi hubungan non-linier melalui fungsi link.
Aplikasi Praktis:
Regresi Logistik: Digunakan dalam analisis risiko penyakit (misalnya, memprediksi risiko diabetes berdasarkan usia dan BMI).
Regresi Poisson: Digunakan dalam analisis kejadian seperti jumlah kecelakaan lalu lintas per hari di suatu kota.
Contoh soal Regresi Logistik:
Kita akan memuat dataset iris dan membuat variabel respons biner: 1 untuk Iris setosa dan 0 untuk spesies lainnya (versicolor atau virginica).
# Memuat dataset iris
data(iris)
# Membuat variabel biner: setosa (1) vs lainnya (0)
iris$is_setosa <- ifelse(iris$Species == "setosa", 1, 0)
# Memilih variabel yang akan digunakan
data <- iris[, c("is_setosa", "Sepal.Length")]
# Melihat beberapa baris data
head(data)
## is_setosa Sepal.Length
## 1 1 5.1
## 2 1 4.9
## 3 1 4.7
## 4 1 4.6
## 5 1 5.0
## 6 1 5.4
Model Regresi Logistik
Kita akan memodelkan probabilitas bunga termasuk Iris setosa berdasarkan panjang sepal (Sepal.Length). Model regresi logistik didefinisikan sebagai:
\[ \log\left( \frac{\pi_i}{1 - \pi_i} \right) = \beta_0 + \beta_1 \cdot \text{Sepal.Length}_i \]
Di mana \(\pi_i\) adalah probabilitas bunga ke-\(i\) adalah Iris setosa.
Perhitungan Model
# Membuat model regresi logistik
model <- glm(is_setosa ~ Sepal.Length, data = data, family = binomial)
# Ringkasan model
summary(model)
##
## Call:
## glm(formula = is_setosa ~ Sepal.Length, family = binomial, data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 27.8285 4.8276 5.765 8.19e-09 ***
## Sepal.Length -5.1757 0.8934 -5.793 6.90e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 190.954 on 149 degrees of freedom
## Residual deviance: 71.836 on 148 degrees of freedom
## AIC: 75.836
##
## Number of Fisher Scoring iterations: 7
exp(coef(model))
## (Intercept) Sepal.Length
## 1.218353e+12 5.652270e-03
Kesimpulan:
Model GLM:
Nilai intercept sebesar 27,83 artinya log-odds bahwa sebuah bunga adalah setosa ketika Sepal Length sama dengan nol adalah sebesar 27,83. Dengan kata lain, ketika panjang kelopak adalah 0 cm, log-odds untuk menjadi setosa cukup tinggi, yang menunjukkan probabilitas yang cukup besar.
Nilai koefisien Sepal Length sebesar -5,18 artinya, setiap kenaikan satu unit (1 cm) pada Sepal Length, maka log-odds bahwa bunga tersebut adalah setosa menurun sebesar 5,18. Ini menunjukkan bahwa panjang kelopak yang semakin besar, secara signifikan menurunkan log-odds menjadi setosa.
Signifikansi Prediktor: Koefisien panjang sepal (Sepal.Length) signifikan secara statistik (p-value (6,90e-09) < 0.05), menunjukkan bahwa panjang sepal berpengaruh terhadap probabilitas bunga apakah termasuk Iris setosa atau tidak.
Eksponen koefisien model:
Berdasarkan perhitungan eksponen di atas, dapat diketahui bahwa semakin panjang sepal length suatu bunga, maka odds atau peluang bunga tersebut merupakan setosa semakin berkurang.
Visualisasi
# Membuat prediksi probabilitas
data$pred_prob <- predict(model, type = "response")
# Plot
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
ggplot(data, aes(x = Sepal.Length, y = pred_prob)) +
geom_point(aes(color = factor(is_setosa)), size = 3, alpha = 0.6) +
geom_line(color = "blue", size = 1) +
labs(title = "Probabilitas Prediksi Spesies Setosa berdasarkan Sepal.Length",
x = "Sepal Length (cm)",
y = "Predicted Probability of Setosa",
color = "Actual Class") +
scale_color_manual(values = c("red", "green"), labels = c("Not Setosa", "Setosa")) +
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.
# Membuat prediksi kelas berdasarkan threshold 0.5
data$pred_class <- ifelse(data$pred_prob > 0.5, 1, 0)
# Membuat confusion matrix
conf_matrix <- table(Predicted = data$pred_class, Actual = data$is_setosa)
conf_matrix
## Actual
## Predicted 0 1
## 0 94 10
## 1 6 40
# Menghitung akurasi
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat("Akurasi:", accuracy, "\n")
## Akurasi: 0.8933333
Kesimpulan:
Visualisasi:
Titik biru menunjukkan bunga yang termasuk Iris setosa, sedangkan titik merah menunjukkan bunga yang bukan Iris setosa.
Garis hitam adalah kurva probabilitas yang diprediksi, menunjukkan bahwa probabilitas menjadi Iris setosa meningkat seiring dengan penurunan panjang sepal.
Evaluasi model:
Berdasarkan perhitungan confusion matrix di atas, dapat diketahui bahwa model tersebut dapat memprediksi dengan cukup akurat dengan jumlah pengelompokkan benar sebanyak 134 dari 150 sampel.
Contoh soal Regresi Poisson:
Sebuah studi ingin mengetahui pengaruh gender dan tingkat kemampuan terhadap jumlah presensi siswa. Berikut data siswa-siswa tersebut:
Gunakan regresi Poisson untuk menganalisis pengaruh gender dan tingkat kemampuan terhadap jumlah presensi.
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
head(quine)
## Eth Sex Age Lrn Days
## 1 A M F0 SL 2
## 2 A M F0 SL 11
## 3 A M F0 SL 14
## 4 A M F0 AL 5
## 5 A M F0 AL 5
## 6 A M F0 AL 13
poisson_model <- glm(Days ~ Sex + Lrn, family = poisson, data = quine)
summary(poisson_model)
##
## Call:
## glm(formula = Days ~ Sex + Lrn, family = poisson, data = quine)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.66232 0.03610 73.752 < 2e-16 ***
## SexM 0.18289 0.04130 4.429 9.49e-06 ***
## LrnSL 0.11777 0.04147 2.840 0.00451 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2073.5 on 145 degrees of freedom
## Residual deviance: 2049.2 on 143 degrees of freedom
## AIC: 2643.7
##
## Number of Fisher Scoring iterations: 5
exp(coef(poisson_model))
## (Intercept) SexM LrnSL
## 14.329477 1.200681 1.124990
quine$predicted <- predict(poisson_model, type = "response")
library(ggplot2)
ggplot(quine, aes(x = Lrn, y=Days, color = Sex)) +
geom_jitter(width = 0.2, alpha = 0.6) +
geom_point(aes(y =predicted), shape = 18, size = 3, color = "black") +
facet_wrap(~Sex) +
labs(title = "Prediksi Jumlah Patahan Benang berdasarkan Wol dan Ketegangan",
x="Tingkat Ketegangan",
y="Jumlah Patahan Benang",
color = "Jenis Wol") +
theme_minimal()
plot(poisson_model$residuals,
main = "Residual Plot",
ylab = "Residual",
xlab = "Index",
pch = 19,
col = "blue")
abline(h=0, col = "red", lty = 2)
Kesimpulan
Model regresi Poisson:
Nilai intercept sebesar 2.66232 artinya log dari rata-rata jumlah hari absen untuk siswa perempuan (Sex = F) dengan kemampuan belajar normal (Lrn = AL) adalah 2.66232 ketika semua prediktor bernilai nol.
Nilai koefisien SexM sebesar 0.18289 artinya, untuk siswa laki-laki dibandingkan perempuan, log dari rata-rata jumlah hari absen meningkat sebesar 0.18289, dengan asumsi kemampuan belajar tetap.
Nilai koefisien LrnSL sebesar 0.11777 artinya, untuk siswa dengan kemampuan belajar lambat (Lrn = SL) dibandingkan siswa dengan kemampuan normal, log dari rata-rata jumlah hari absen meningkat sebesar 0.11777, dengan asumsi jenis kelamin tetap.
p-value untuk SexM adalah 9.49e-06 (< 0.05), menunjukkan bahwa jenis kelamin (Sex) adalah prediktor yang signifikan terhadap jumlah hari absen. Siswa laki-laki cenderung memiliki hari absen yang lebih banyak dibandingkan siswa perempuan.
p-value untuk LrnSL adalah 0.00451 (< 0.05), menunjukkan bahwa kemampuan belajar (Lrn) juga merupakan prediktor yang signifikan. Siswa dengan kemampuan belajar lambat cenderung memiliki hari absen yang lebih banyak dibandingkan siswa dengan kemampuan belajar normal.
Residual Deviance: 2049.2 dengan 143 derajat kebebasan, menunjukkan bahwa model cukup baik dalam menjelaskan variabilitas data, meskipun masih ada deviasi yang tidak dijelaskan.
AIC: 2643.7, yang dapat digunakan untuk membandingkan model ini dengan model lain. Nilai AIC yang lebih kecil menunjukkan model yang lebih baik.
Eksponensial koefisien model:
Nilai untuk intercept adalah \(e^{2.66232} \approx 14.329477\). Ini artinya rata-rata jumlah hari absen untuk siswa perempuan dengan kemampuan belajar normal adalah sekitar 14.33 hari.
Nilai untuk SexM adalah \(e^{0.18289} \approx 1.200681\). Ini artinya, rata-rata jumlah hari absen untuk siswa laki-laki meningkat sebesar 1.200681 kali lipat (atau sekitar 20% lebih tinggi) dibandingkan siswa perempuan, dengan kemampuan belajar yang sama.
Nilai untuk LrnSL adalah \(e^{0.11777} \approx 1.124990\). Ini artinya, rata-rata jumlah hari absen untuk siswa dengan kemampuan belajar lambat meningkat sebesar 1.124990 kali lipat (atau sekitar 12.5% lebih tinggi) dibandingkan siswa dengan kemampuan normal, dengan jenis kelamin yang sama.
Generalized Linear Model (GLM) adalah kerangka statistik yang digunakan untuk memodelkan hubungan antara variabel respons dan prediktor, terutama untuk data yang tidak mengikuti distribusi normal, seperti data biner (regresi logistik) atau data hitung (regresi Poisson). GLM terdiri dari tiga komponen utama: distribusi acak (dari keluarga eksponensial), komponen sistematik (kombinasi linier prediktor), dan fungsi link yang menghubungkan keduanya. Bagian ini membahas ekspektasi, varians, metode penaksiran, diagnostik, serta estimasi dan inferensi untuk regresi logistik dan Poisson.
Dalam GLM, variabel respons \(Y\) diasumsikan mengikuti distribusi dari keluarga eksponensial, seperti Bernoulli untuk data biner atau Poisson untuk data hitung. Ekspektasi dan varians \(Y\) bergantung pada distribusi yang dipilih.
Ekspektasi: Ekspektasi \(Y\), dilambangkan sebagai \(\mu\), adalah rata-rata dari distribusi \(Y\). Dalam GLM, \(\mu\) dihubungkan dengan prediktor melalui fungsi link: \(g(\mu) = \mathbf{x}^\top \beta\), di mana \(\mathbf{x}\) adalah vektor prediktor dan \(\beta\) adalah parameter.
Varians: Varians \(Y\) biasanya merupakan fungsi dari ekspektasi \(\mu\), dinyatakan sebagai \(V(Y) = \phi v(\mu)\), di mana \(\phi\) adalah parameter dispersi dan \(v(\mu)\) adalah fungsi varians yang spesifik untuk distribusi.
Contoh Distribusi
Regresi Logistik (Bernoulli): \(Y \sim \text{Bernoulli}(\pi)\), dengan \(\mu = \pi\), \(V(Y) = \pi(1-\pi)\), dan \(\phi = 1\).
Regresi Poisson: \(Y \sim \text{Poisson}(\lambda)\), dengan \(\mu = \lambda\), \(V(Y) = \lambda\), dan \(\phi = 1\).
Untuk menghitung ekspektasi dan varians dalam GLM, kita perlu menentukan distribusi \(Y\) dan fungsi link yang sesuai.
Regresi Logistik
Distribusi: \(Y_i \sim \text{Bernoulli}(\pi_i)\).
Ekspektasi: \(\mu_i = \pi_i = P(Y_i = 1 | \mathbf{x}_i) = \frac{\exp(\mathbf{x}_i^\top \beta)}{1 + \exp(\mathbf{x}_i^\top \beta)}\).
Varians: \(V(Y_i) = \pi_i (1 - \pi_i)\).
Fungsi Link: Logit, \(g(\mu_i) = \log\left(\frac{\mu_i}{1 - \mu_i}\right) = \mathbf{x}_i^\top \beta\).
Regresi Poisson
Distribusi: \(Y_i \sim \text{Poisson}(\lambda_i)\).
Ekspektasi: \(\mu_i = \lambda_i = \exp(\mathbf{x}_i^\top \beta)\).
Varians: \(V(Y_i) = \lambda_i = \mu_i\).
Fungsi Link: Log, \(g(\mu_i) = \log(\mu_i) = \mathbf{x}_i^\top \beta\).
Parameter \(\beta\) dalam GLM biasanya ditaksir menggunakan Maximum Likelihood Estimation (MLE). Fungsi likelihood untuk \(n\) observasi adalah:
\[ L(\beta) = \prod_{i=1}^n f(y_i | \mathbf{x}_i, \beta) \]
Log-likelihood:
\[ \ell(\beta) = \sum_{i=1}^n \log f(y_i | \mathbf{x}_i, \beta) \]
Karena log-likelihood sering kali non-linear, metode iteratif seperti Newton-Raphson atau Iteratively Reweighted Least Squares (IRLS) digunakan untuk menemukan \(\hat{\beta}\).
Newton-Raphson
Iterasi Newton-Raphson didefinisikan sebagai:
\[ \beta^{(t+1)} = \beta^{(t)} - H^{-1}(\beta^{(t)}) U(\beta^{(t)}) \]
Di mana:
\(U(\beta) = \frac{\partial \ell(\beta)}{\partial \beta}\) adalah score function.
\(H(\beta) = \frac{\partial^2 \ell(\beta)}{\partial \beta \partial \beta^\top}\) adalah Hessian matrix.
Diagnostik model bertujuan untuk mengevaluasi kecocokan model GLM terhadap data. Beberapa metode diagnostik meliputi:
Statistik Devians
Devians mengukur perbedaan antara model yang ditaksir dan model jenuh (saturated model):
\[ D = 2 \sum \left[ y_i \log\left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i - \hat{\mu}_i) \right] \]
Devians kecil menunjukkan model yang lebih cocok.
Statistik Chi-Square Pearson
Statistik Pearson mengukur perbedaan antara nilai observasi dan prediksi:
\[ X^2 = \sum \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i} \]
Nilai kecil menunjukkan kecocokan yang baik.
Analisis Residual
Residual dihitung sebagai \(r_i = y_i - \hat{\mu}_i\). Plot residual terhadap indeks atau nilai prediksi membantu mendeteksi penyimpangan sistematis.
Estimasi Parameter
Dalam regresi logistik, \(Y_i \sim \text{Bernoulli}(\pi_i)\), dengan:
\[ \pi_i = \frac{\exp(\mathbf{x}_i^\top \beta)}{1 + \exp(\mathbf{x}_i^\top \beta)} \]
Log-likelihood:
\[ \ell(\beta) = \sum_{i=1}^n \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right] \]
Score Function:
\[ U(\beta) = \mathbf{X}^\top (\mathbf{y} - \pi) \]
Hessian Matrix:
\[ H(\beta) = -\mathbf{X}^\top \mathbf{W} \mathbf{X}, \quad \mathbf{W} = \text{diag}(\pi_i (1 - \pi_i)) \]
Iterasi Newton-Raphson:
\[ \beta^{(t+1)} = \beta^{(t)} + (\mathbf{X}^\top \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^\top (\mathbf{y} - \pi^{(t)}) \]
Inferensi
Uji Wald
Uji Wald menguji signifikansi parameter \(\beta_j\):
\[ z = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)} \sim \mathcal{N}(0,1) \]
Uji Likelihood Ratio (LRT)
LRT membandingkan model penuh dan model nol:
\[ \Lambda = -2 \left( \ell(\beta_{\text{nol}}) - \ell(\beta_{\text{penuh}}) \right) \sim \chi^2_{\text{df}} \]
Evaluasi Model
Model dievaluasi menggunakan AIC dan BIC:
AIC: \(-2 \ell(\hat{\beta}) + 2p\)
BIC: \(-2 \ell(\hat{\beta}) + p \log(n)\)
Estimasi Parameter
Dalam regresi Poisson, \(Y_i \sim \text{Poisson}(\lambda_i)\), dengan:
\[ \lambda_i = \exp(\mathbf{x}_i^\top \beta) \]
Log-likelihood:
\[ \ell(\beta) = \sum_{i=1}^n \left[ y_i \log(\lambda_i) - \lambda_i - \log(y_i!) \right] \]
Score Function:
\[ U(\beta) = \mathbf{X}^\top (\mathbf{y} - \lambda) \]
Hessian Matrix:
\[ H(\beta) = -\mathbf{X}^\top \mathbf{W} \mathbf{X}, \quad \mathbf{W} = \text{diag}(\lambda_i) \]
Iterasi IRLS (ekivalen dengan Newton-Raphson):
\[ \beta^{(t+1)} = (\mathbf{X}^\top \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{W}^{(t)} \mathbf{z}^{(t)} \]
Di mana \(\mathbf{z}^{(t)} = \eta^{(t)} + \frac{\mathbf{y} - \lambda^{(t)}}{\lambda^{(t)}}\), dan \(\eta^{(t)} = \mathbf{X} \beta^{(t)}\).
Inferensi
Uji Wald
Sama seperti regresi logistik, uji Wald digunakan untuk \(\beta_j\):
\[ z = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)} \sim \mathcal{N}(0,1) \]
Uji Likelihood Ratio (LRT)
LRT juga berlaku untuk regresi Poisson:
\[ \Lambda = -2 \left( \ell(\beta_{\text{nol}}) - \ell(\beta_{\text{penuh}}) \right) \sim \chi^2_{\text{df}} \]
Evaluasi Model
Sama seperti regresi logistik, AIC dan BIC digunakan untuk mengevaluasi model.
Contoh soal metode estimasi dan inferensi regresi logistik:
Seorang peneliti ingin memodelkan probabilitas mahasiswa diterima di program beasiswa (diterima = 1, tidak diterima = 0) berdasarkan nilai rata-rata akademik mereka (\(x\), dalam skala 0-100). Data disimulasikan untuk 200 mahasiswa dengan log-odds keberhasilan \(-3 + 0.05x\). Gunakan metode Newton-Raphson untuk estimasi parameter, lalu lakukan uji Wald dan uji Likelihood Ratio (Chi-Square) untuk inferensi.
Penyelesaian
Langkah 1: Simulasi Data
set.seed(123)
x <- rnorm(200, mean=75, sd=10) # Nilai akademik (mean=75, sd=10)
X <- cbind(1, x)
beta_true <- c(-3, 0.05)
eta <- X %*% beta_true
p <- 1 / (1 + exp(-eta))
y <- rbinom(200, 1, p)
Langkah 2: Estimasi dengan Newton-Raphson
Kita gunakan metode Newton-Raphson untuk mengestimasi parameter \(\beta\).
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) {
cat("Konvergen pada iterasi ke-", i, "\n")
break
}
beta <- beta_new
}
## Konvergen pada iterasi ke- 5
cat("Estimasi parameter:\n")
## Estimasi parameter:
print(beta)
## [,1]
## -1.30630586
## x 0.02735003
Langkah 3: Inferensi
Uji Wald
Kita lakukan uji Wald untuk menguji signifikansi prediktor \(x\).
Hipotesis Uji:
\(H_0: \beta_1 = 0\) (nilai akademik tidak berpengaruh terhadap probabilitas diterima).
\(H_1: \beta_1 \neq 0\) (nilai akademik berpengaruh terhadap probabilitas diterima).
model <- glm(y ~ x, family = binomial)
summary_model <- summary(model)
coef_val <- coef(model)[2]
se_val <- summary_model$coefficients[2, 2]
wald_z <- coef_val / se_val
p_value <- 2 * (1 - pnorm(abs(wald_z)))
cat("Statistik Wald: Z =", wald_z, "\np-value =", p_value, "\n")
## Statistik Wald: Z = 1.642775
## p-value = 0.1004294
Keputusan:
Karena p-value > 0.05, maka terima \(H_0\). Nilai akademik tidak berpengaruh signifikan terhadap probabilitas diterima di program beasiswa.
Uji Likelihood Ratio (Chi-Square)
Kita bandingkan model penuh dengan model nol menggunakan uji LRT.
Hipotesis Uji:
\(H_0\): Model nol cukup (tidak ada prediktor yang signifikan).
\(H_1\): Model penuh lebih baik (setidaknya satu prediktor signifikan).
model_null <- glm(y ~ 1, family = binomial)
lrt <- anova(model_null, model, test = "Chisq")
cat("Uji Likelihood Ratio:\n")
## Uji Likelihood Ratio:
print(lrt)
## Analysis of Deviance Table
##
## Model 1: y ~ 1
## Model 2: y ~ x
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 199 252.23
## 2 198 249.45 1 2.7845 0.09518 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Keputusan:
Karena p-value > 0.05, maka terima \(H_0\). Model nol lebih baik dibanding model penuh dengan prediktor nilai akademik.
Kesimpulan
Estimasi parameter menggunakan Newton-Raphson menghasilkan nilai yang mendekati nilai sebenarnya (\(\beta_0 = -3\), \(\beta_1 = 0.05\)). Uji Wald menunjukkan p-value (> 0.05), sehingga kita menerima \(H_0\) dan menyimpulkan bahwa nilai akademik (\(x\)) tidak berpengaruh signifikan terhadap probabilitas diterima di program beasiswa. Uji LRT juga menghasilkan p-value (> 0.05), mengindikasikan bahwa model nol lebih baik dibandingkan dengan model prediktor.
Contoh Soal metode estimasi dan inferensi regresi poisson:
Seorang peneliti memodelkan jumlah pesanan harian di sebuah restoran berdasarkan suhu harian (\(x\), dalam derajat Celsius). Data disimulasikan untuk 180 hari dengan \(\lambda = \exp(2.2 + 0.03x)\). Gunakan metode Iteratively Reweighted Least Squares (IRLS) untuk estimasi parameter, lalu lakukan uji Wald dan uji Likelihood Ratio (Chi-Square) untuk inferensi.
Penyelesaian
Langkah 1: Simulasi Data
set.seed(456)
x <- rnorm(180, mean=25, sd=5) # Suhu harian (mean=25, sd=5)
X <- cbind(1, x)
beta_true <- c(2.2, 0.03)
eta <- X %*% beta_true
lambda <- exp(eta)
y <- rpois(180, lambda)
Langkah 2: Estimasi dengan IRLS
Kita gunakan metode IRLS untuk mengestimasi parameter \(\beta\).
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- 26
cat("Hasil estimasi akhir:\n")
## Hasil estimasi akhir:
print(beta)
## [,1]
## 2.14663531
## x 0.03266819
Langkah 3: Inferensi
Uji Wald
Kita lakukan uji Wald untuk menguji signifikansi prediktor \(x\).
Hipotesis Uji:
\(H_0: \beta_1 = 0\) (suhu harian tidak berpengaruh terhadap jumlah pesanan).
\(H_1: \beta_1 \neq 0\) (suhu harian berpengaruh terhadap jumlah pesanan).
data <- data.frame(x = x, y = y)
model <- glm(y ~ x, family = poisson, data = data)
summary_model <- summary(model)
coef_val <- coef(model)[2]
se_val <- summary_model$coefficients[2, 2]
wald_z <- coef_val / se_val
p_value <- 2 * (1 - pnorm(abs(wald_z)))
cat("Statistik Wald: Z =", wald_z, "\np-value =", p_value, "\n")
## Statistik Wald: Z = 9.697236
## p-value = 0
Keputusan:
Karena p-value < 0.05, maka tolak \(H_0\). Suhu harian berpengaruh signifikan terhadap jumlah pesanan harian di restoran.
Uji Likelihood Ratio (Chi-Square)
Kita bandingkan model penuh dengan model nol menggunakan uji LRT.
Hipotesis Uji:
\(H_0\): Model nol cukup (tidak ada prediktor yang signifikan).
\(H_1\): Model penuh lebih baik (setidaknya satu prediktor signifikan).
model_null <- glm(y ~ 1, family = poisson, data = data)
lrt <- anova(model_null, model, test = "Chisq")
cat("Uji Likelihood Ratio:\n")
## Uji Likelihood Ratio:
print(lrt)
## Analysis of Deviance Table
##
## Model 1: y ~ 1
## Model 2: y ~ x
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 179 264.56
## 2 178 170.06 1 94.506 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Keputusan:
Karena p-value < 0.05, maka tolak \(H_0\). Model penuh dengan prediktor suhu harian lebih baik dibandingkan model nol.
Kesimpulan
Estimasi parameter menggunakan IRLS menghasilkan nilai yang mendekati nilai sebenarnya (\(\beta_0 = 2.2\), \(\beta_1 = 0.03\)). Uji Wald menunjukkan p-value kecil (< 0.05), sehingga kita menolak \(H_0\) dan menyimpulkan bahwa suhu harian (\(x\)) berpengaruh signifikan terhadap jumlah pesanan harian di restoran. Uji LRT juga menghasilkan p-value kecil, mengindikasikan bahwa model dengan prediktor lebih baik dibandingkan model nol, sehingga prediktor \(x\) berkontribusi signifikan pada model.
Deskripsi data
Data berikut didapat dari website archive.ics.uc.edu. Data tersebut berisi tentang kampanye pemasaran langsung (panggilan telepon) dari sebuah lembaga perbankan Portugis. Tujuan klasifikasi adalah untuk memprediksi apakah klien akan berlangganan deposito berjangka atau tidak (variabel y) jika terdapat beberapa variabel prediktor yang dikontrol.
# Load the data
library(openxlsx)
## Warning: package 'openxlsx' was built under R version 4.3.3
data <- read.xlsx(file.choose())
# Keep only numeric variables and the target variable 'y'
numeric_vars <- c("age", "balance", "day", "duration", "campaign", "pdays", "previous", "y")
data <- data[, numeric_vars]
# Convert 'y' to binary (1 for yes, 0 for no)
data$y <- ifelse(data$y == "yes", 1, 0)
head(data)
## age balance day duration campaign pdays previous y
## 1 30 1787 19 79 1 -1 0 0
## 2 33 4789 11 220 1 339 4 0
## 3 35 1350 16 185 1 330 1 0
## 4 30 1476 3 199 4 -1 0 0
## 5 59 0 5 226 1 -1 0 0
## 6 35 747 23 141 2 176 3 0
# Check the structure of the data
str(data)
## 'data.frame': 4521 obs. of 8 variables:
## $ age : num 30 33 35 30 59 35 36 39 41 43 ...
## $ balance : num 1787 4789 1350 1476 0 ...
## $ day : num 19 11 16 3 5 23 14 6 14 17 ...
## $ duration: num 79 220 185 199 226 141 341 151 57 313 ...
## $ campaign: num 1 1 1 4 1 2 1 2 2 1 ...
## $ pdays : num -1 339 330 -1 -1 176 330 -1 -1 147 ...
## $ previous: num 0 4 1 0 0 3 2 0 0 2 ...
## $ y : num 0 0 0 0 0 0 0 0 0 0 ...
# Fit the logistic regression model
model <- glm(y ~ age + balance + day + duration + campaign + pdays + previous,
data = data,
family = binomial(link = "logit"))
# Summary of the model (includes Wald test for coefficients)
summary(model)
##
## Call:
## glm(formula = y ~ age + balance + day + duration + campaign +
## pdays + previous, family = binomial(link = "logit"), data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.927e+00 2.498e-01 -15.724 < 2e-16 ***
## age 1.357e-02 4.674e-03 2.903 0.003694 **
## balance 2.478e-05 1.484e-05 1.670 0.095003 .
## day 5.019e-03 6.561e-03 0.765 0.444244
## duration 3.644e-03 1.758e-04 20.728 < 2e-16 ***
## campaign -9.074e-02 2.600e-02 -3.490 0.000483 ***
## pdays 1.882e-03 5.066e-04 3.716 0.000203 ***
## previous 1.044e-01 2.626e-02 3.976 7.01e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 3231.0 on 4520 degrees of freedom
## Residual deviance: 2613.4 on 4513 degrees of freedom
## AIC: 2629.4
##
## Number of Fisher Scoring iterations: 6
Kesimpulan:
Intersep -3,93 menunjukkan log-odds probabilitas berlangganan deposito jika semua variabel prediktor bernilai nol
Koefisien 1,357e-02 berarti bahwa untuk setiap kenaikan 1 umur pelanggan, log-odds probabilitas berlangganan deposito meningkat sebesar 1,357e-02, dengan asumsi variabel lain konstan.
Koefisien 2.478e-05 berarti bahwa untuk setiap kenaikan 1 euro pada saldo pelanggan, log-odds probabilitas berlangganan deposito meningkat sebesar 0.00002478, dengan asumsi variabel lain konstan.
Koefisien 5.019e-03 berarti bahwa untuk setiap kenaikan 1 bulan ketika klien terakhir kali dihubungi, log-odds probabilitas berlangganan deposito meningkat sebesar 5.019e-03, dengan asumsi variabel lain konstan.
Koefisien 3.644e-03 berarti bahwa untuk setiap kenaikan 1 lamanya durasi kontak dengan klien, log-odds probabilitas berlangganan deposito meningkat sebesar 3.644e-03, dengan asumsi variabel lain konstan.
Koefisien -9.074e-02 berarti bahwa untuk setiap kenaikan 1 kontak dengan klien ketika kampanye, log-odds probabilitas berlangganan deposito menurun sebesar -9.074e-02, dengan asumsi variabel lain konstan.
Koefisien 1.882e-03 berarti bahwa untuk setiap kenaikan 1 hari sejak klien terakhir dihubungi, log-odds probabilitas berlangganan deposito meningkat sebesar 1.882e-03, dengan asumsi variabel lain konstan.
Koefisien 1.044e-01 berarti bahwa untuk setiap kenaikan 1 kontak kepada klien sebelum kampanye, log-odds probabilitas berlangganan deposito meningkat sebesar 1.044e-01, dengan asumsi variabel lain konstan.
Koefisien positif pada age, balance, day, duration, pdays, dan previous menunjukkan bahwa peningkatan nilai variabel ini juga akan meningkatkan log-odds kemungkinan berlangganan. Sebaliknya, koefisien negatif pada campaign menunjukkan bahwa peningkatan jumlah kontak kampanye justru menurunkan log-odds kemungkinan berlangganan.
AIC sebesar 2629.4 dapat digunakan untuk membandingkan model ini dengan model lain (nilai AIC yang lebih kecil menunjukkan model yang lebih baik).
Hipotesis Uji:
\(H_0\): Model nol cukup (tidak ada prediktor yang signifikan).
\(H_1\): Model penuh lebih baik (setidaknya satu prediktor signifikan).
# Fit the null model (intercept only)
null_model <- glm(y ~ 1, data = data, family = binomial(link = "logit"))
# Compute the Likelihood Ratio Test statistic
lr_stat <- 2 * (logLik(model) - logLik(null_model))
df <- length(coef(model)) - 1 # Degrees of freedom
p_value_lr <- 1 - pchisq(lr_stat, df)
# Print results
cat("Likelihood Ratio Test:\n")
## Likelihood Ratio Test:
cat("LR Statistic:", lr_stat, "\n")
## LR Statistic: 617.6131
cat("P-value:", p_value_lr, "\n")
## P-value: 0
Kesimpulan
P-value = 0 (sangat kecil, < 0.05) menunjukkan bahwa kita menolak hipotesis nol. Ini berarti model penuh dengan prediktor (age, balance, day, duration, campaign, pdays, previous) secara signifikan lebih baik dalam menjelaskan variabilitas probabilitas seorang pelanggan berlangganan depsito dibandingkan model nol. Dengan kata lain, setidaknya satu dari prediktor memiliki pengaruh signifikan terhadap probabilitas seorang pelanggan berlangganan depsito.
Hipotesis Uji:
\(H_0: \beta_1 = 0\) (variabel prediktor tidak berpengaruh terhadap probabilitas berlangganan deposito).
\(H_1: \beta_1 \neq 0\) (variabel prediktor berpengaruh terhadap probabilitas berlangganan deposito).
# Extract Wald test p-values
wald_pvalues <- summary(model)$coefficients[, 4]
cat("Wald Test P-values for Coefficients:\n")
## Wald Test P-values for Coefficients:
print(wald_pvalues)
## (Intercept) age balance day duration campaign
## 1.040022e-55 3.693962e-03 9.500296e-02 4.442444e-01 1.950557e-95 4.833191e-04
## pdays previous
## 2.025523e-04 7.011295e-05
Kriteria uji
Tolak H0 jika p-value < 0,05
Kesimpulan
Variabel age, duration, campaign, pdays, dan previous signifikan (p-value < 0.05), sehingga masing-masing variabel ini memiliki pengaruh yang signifikan terhadap probabilitas seorang pelanggan berlangganan deposito (setelah mengontrol variabel lain).
Variabel balance (p-value = 0.095003) dan day (p-value = 0.444244) tidak signifikan (p-value > 0.05), sehingga tidak ada bukti yang cukup bahwa kedua variabel ini memengaruhi probabilitas seorang pelanggan berlangganan deposito.
Regresi logistik merupakan model statistik yang digunakan untuk memodelkan hubungan antara satu variabel respon kategorik biner (seperti “ya/tidak”, “lulus/gagal”) dengan satu atau lebih variabel prediktor. Dalam praktiknya, variabel prediktor dapat berasal dari berbagai jenis skala pengukuran:
Nominal: Kategori yang tidak memiliki urutan (misalnya: jenis kelamin, tempat tinggal).
Ordinal: Kategori yang memiliki urutan, tapi tidak memiliki jarak tetap (misalnya: tingkat pendidikan).
Rasio: Variabel kuantitatif dengan nol absolut dan jarak tetap (misalnya: penghasilan, usia, jumlah jam belajar).
Ketiga jenis prediktor tersebut dapat dimasukkan dalam model regresi logistik, dengan perlakuan yang berbeda tergantung jenis skalanya:
Prediktor nominal dikodekan menggunakan dummy variabel agar dapat diproses dalam model.
Prediktor ordinal bisa diperlakukan sebagai nominal (dummy) atau sebagai numeric rank tergantung asumsi yang digunakan.
Prediktor rasio bisa langsung digunakan dalam bentuk aslinya tanpa transformasi.
Dikutip dari Agresti (2010), perlakuan terhadap variabel ordinal dalam regresi logistik harus memperhatikan konteks dan tujuan analisis:
Jika urutan kategori penting dan diasumsikan memiliki efek linear, maka perlakuan sebagai numeric rank lebih tepat.
Namun, jika urutan tidak diasumsikan memberi jarak yang sama, lebih baik gunakan pendekatan dummy.
Interpretasi koefisien dari regresi logistik dilakukan dalam skala log-odds. Namun, untuk memudahkan interpretasi, koefisien tersebut biasanya dikonversi ke odds ratio dengan cara mengaplikasikan fungsi eksponensial pada koefisien model.
Simulasi Data
library(tibble)
set.seed(2025)
n <- 600
jenis_kelamin <- sample(c("Laki-laki", "Perempuan"), n, replace = TRUE)
tingkat_pendidikan <- sample(c("SMA", "D3", "S1", "S2"), n, replace = TRUE, prob = c(0.3, 0.3, 0.25, 0.15))
jumlah_jam_belajar <- round(rnorm(n, mean = 12, sd = 4), 1)
jumlah_jam_belajar[jumlah_jam_belajar < 0] <- 0 # tidak boleh negatif
# Buat logit dan probabilitas
logit_p <- -3 + 0.6 * (jenis_kelamin == "Perempuan") + 0.7 * as.numeric(factor(tingkat_pendidikan, ordered = TRUE)) + 0.15 * jumlah_jam_belajar
p <- 1 / (1 + exp(-logit_p))
# Simulasi hasil kelulusan
lulus <- rbinom(n, 1, p)
data_ujian <- tibble(lulus, jenis_kelamin, tingkat_pendidikan, jumlah_jam_belajar)
head(data_ujian)
## # A tibble: 6 × 4
## lulus jenis_kelamin tingkat_pendidikan jumlah_jam_belajar
## <int> <chr> <chr> <dbl>
## 1 1 Laki-laki SMA 8.6
## 2 0 Perempuan S1 5.9
## 3 1 Perempuan D3 16.1
## 4 1 Perempuan D3 21.2
## 5 1 Laki-laki D3 14.8
## 6 1 Laki-laki D3 12.3
Eksplorasi Data
data_ujian %>%
dplyr::group_by(lulus) %>%
dplyr::summarise(
Jumlah = dplyr::n(),
jumlahjambelajar = mean(jumlah_jam_belajar)
)
## # A tibble: 2 × 3
## lulus Jumlah jumlahjambelajar
## <int> <int> <dbl>
## 1 0 229 10.9
## 2 1 371 12.6
Berdasarkan perhitungan eksplorasi data tersebut, dapat diketahui jumlah lulus ada sebanyak 371 dengan rata-rata jam belajar sebesar 12,598 jam. Sedangkann jumlah yang tidak lulus ada sebanyak 229 dengan rata-rata jam belajar sebesar 10,851 jam.
Perlakuan Variabel Ordinal
library(dplyr)
data_ujian_nominal <- data_ujian %>%
mutate(
education = factor(tingkat_pendidikan, levels = c("SMA", "D3", "S1", "S2"))
)
model_nominal <- glm(lulus ~ jenis_kelamin + tingkat_pendidikan + jumlah_jam_belajar,
family = binomial, data = data_ujian_nominal)
summary(model_nominal)
##
## Call:
## glm(formula = lulus ~ jenis_kelamin + tingkat_pendidikan + jumlah_jam_belajar,
## family = binomial, data = data_ujian_nominal)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.12723 0.35437 -6.003 1.94e-09 ***
## jenis_kelaminPerempuan 0.51680 0.18891 2.736 0.00622 **
## tingkat_pendidikanS1 0.75097 0.23640 3.177 0.00149 **
## tingkat_pendidikanS2 1.40712 0.27969 5.031 4.88e-07 ***
## tingkat_pendidikanSMA 2.24781 0.27183 8.269 < 2e-16 ***
## jumlah_jam_belajar 0.12107 0.02493 4.856 1.20e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 797.85 on 599 degrees of freedom
## Residual deviance: 667.39 on 594 degrees of freedom
## AIC: 679.39
##
## Number of Fisher Scoring iterations: 4
Keterangan: Model menggunakan education sebagai variabel dummy, baseline “D3”. Setiap koefisien membandingkan kategori terhadap baseline.
Interpretasi: Intercept: -2,127 - Ini adalah log-odds dasar untuk individu berjenis kelamin laki-laki, berpendidikan D3, dan memiliki jumlah jam belajar sebesar 0. - Signifikan (p = 1,94e-09), artinya baseline ini berbeda secara signifikan dari probabilitas 50%. - Odds Ratio = 0,114 menunjukkan bahwa peluang sukses lebih besar daripada baseline.
Jenis_kelaminPerempuan: 0.51680 - Individu dengan jenis kelamin perempuan memiliki log-odds sukses 0,5 lebih tinggi dibandingkan individu dengan jenis kelamin laki-laki. - Signifikan di level 5% dengan nilai p-value = 0.00622, artinya ada indikasi perempuan memiliki peluang sukses yang lebih tinggi dibandingkan laki-laki. - Odds Ratio = 1,677 menunjukkan bahwa peluang sukses perempuan sekitar 66.7% lebih besar dibanding laki-laki.
Tingkat pendidikan S1: 0.75097 - Individu dengan pendidikan S1 memiliki log-odds sukses 0,75 lebih tinggi dibandingkan individu dengan pendidikan D3. - Signifikan di level 5% dengan nilai p-value = 0.00149, artinya ada indikasi individu dengan pendidikan S1 memiliki peluang sukses yang lebih tinggi dibandingkan individu dengan pendidikan D3. - Odds Ratio = 2,117 menunjukkan bahwa peluang sukses individu dengan pendidikan S1 sekitar 2 kali lebih besar dibanding individu dengan pendidikan D3.
Tingkat pendidikan S2: 1.40712 - Individu dengan pendidikan S2 memiliki log-odds sukses 1,41 lebih tinggi dibandingkan individu dengan pendidikan D3. - Signifikan di level 5% dengan nilai p-value = 4.88e-07, artinya ada indikasi individu dengan pendidikan S2 memiliki peluang sukses yang lebih tinggi dibandingkan individu dengan pendidikan D3. - Odds Ratio = 4,096 menunjukkan bahwa peluang sukses individu dengan pendidikan S1 sekitar 4 kali lebih besar dibanding individu dengan pendidikan D3.
Tingkat pendidikan SMA: 2.24781 - Individu dengan pendidikan SMA memiliki log-odds sukses 2,25 lebih tinggi dibandingkan individu dengan pendidikan D3. - Signifikan di level 5% dengan nilai p-value = 2e-16, artinya ada indikasi individu dengan pendidikan SMA memiliki peluang sukses yang lebih tinggi dibandingkan individu dengan pendidikan D3. - Odds Ratio = 9,49 menunjukkan bahwa peluang sukses individu dengan pendidikan SMA sekitar 9 kali lebih besar dibanding individu dengan pendidikan D3.
Jumlah jam belajar: 0.12107 - Setiap kenaikan 1 jam jumlah jam belajar, artinya log-odds lulus meningkat sebesar 0.121. - Signifikan di level 5% dengan nilai p-value = 1.20e-06, artinya jumlah jam belajar berhubungan positif dengan peluang lulus. - Odds Ratio = 1,129 menunjukkan bahwa setiap tambahan 1 jam jumlah jam belajar, maka peluang sukses meningkst sekitar 12,9%
Interpretasi Goodness-of-Fit - Null deviance (797.85): Deviance model tanpa prediktor. - Residual deviance (667.39): Deviance model dengan prediktor. - Penurunan dari null deviance ke residual deviance menunjukkan bahwa model membawa informasi yang cukup baik. - AIC (679.39) digunakan untuk membandingkan model, semakin kecil AIC, semakin baik model dalam menyeimbangkan goodness-of-fit dan kompleksitas.
Signifikansi Model - Variabel jumlah jam belajar dan tingkat pendidikan (S1, S2, SMA) signifikan meningkatkan peluang lulus. - Variabel gender juga signifikan pada taraf 5%.
Kesimpulan Praktis - Jumlah jam belajar dan tingkat pendidikan yang lebih tinggi merupakan prediktor kuat untuk peluang lulus. - Gender juga berpotensi berpengaruh terhadap peluang lulus. - Model cukup baik dalam memprediksi dibandingkan model null.
data_ujian_numeric <- data_ujian %>%
mutate(
education_numeric = case_when(
tingkat_pendidikan == "SMA" ~ 1,
tingkat_pendidikan == "D3" ~ 2,
tingkat_pendidikan == "S1" ~ 3,
tingkat_pendidikan == "S2" ~ 4
)
)
model_numeric <- glm(lulus ~ jenis_kelamin + education_numeric + jumlah_jam_belajar,
family = binomial, data = data_ujian_numeric)
summary(model_numeric)
##
## Call:
## glm(formula = lulus ~ jenis_kelamin + education_numeric + jumlah_jam_belajar,
## family = binomial, data = data_ujian_numeric)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.79089 0.35380 -2.235 0.02539 *
## jenis_kelaminPerempuan 0.61553 0.17512 3.515 0.00044 ***
## education_numeric -0.17250 0.08388 -2.057 0.03973 *
## jumlah_jam_belajar 0.11604 0.02289 5.069 4e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 797.85 on 599 degrees of freedom
## Residual deviance: 752.37 on 596 degrees of freedom
## AIC: 760.37
##
## Number of Fisher Scoring iterations: 4
Interpretasi: Intercept: -0.79089 - Ini adalah log-odds dasar untuk individu berjenis kelamin laki-laki, berpendidikan D3, dan memiliki jumlah jam belajar sebesar 0. - Signifikan (p = 0.02539), artinya baseline ini berbeda secara signifikan dari probabilitas 50%. - Odds Ratio = 0,453 menunjukkan bahwa peluang sukses lebih besar daripada baseline.
Jenis_kelaminPerempuan: 0.61553 - Individu dengan jenis kelamin perempuan memiliki log-odds sukses 0,6 lebih tinggi dibandingkan individu dengan jenis kelamin laki-laki. - Signifikan di level 5% dengan nilai p-value = 0.00044, artinya ada indikasi perempuan memiliki peluang sukses yang lebih tinggi dibandingkan laki-laki. - Odds Ratio = 1,85 menunjukkan bahwa peluang sukses perempuan sekitar 85% lebih besar dibanding laki-laki.
Education numeric: -0.17250 - Setiap kenaikan satu tingkat pendidikan (HighSchool-> Bachelor-> Master-> PhD) menurunkan log-odds lulus sebesar 0.172. - p = 0.03973 (sangat signifikan), menunjukkan bahwa semakin tinggi tingkat pendidikan, justru semakin kecil peluang lulus. - Odds Ratio = 0,842 artinya setiap kenaikan satu tingkat pendidikan menurunkan peluang lulus sekitar 84,2%.
Jumlah jam belajar: 0.11604 - Setiap kenaikan 1 jam jumlah jam belajar, artinya log-odds lulus meningkat sebesar 0.116. - Signifikan di level 5% dengan nilai p-value = 4e-07, artinya jumlah jam belajar berhubungan positif dengan peluang lulus. - Odds Ratio = 1,123 menunjukkan bahwa setiap tambahan 1 jam jumlah jam belajar, maka peluang sukses meningkst sekitar 12,3%.
Interpretasi Goodness-of-Fit - Null deviance (797.85): Deviance model tanpa prediktor. - Residual deviance (752.37): Deviance model dengan prediktor. - Penurunan dari null deviance ke residual deviance menunjukkan bahwa model membawa informasi yang cukup baik. - AIC (760.37) digunakan untuk membandingkan model, semakin kecil AIC, semakin baik model dalam menyeimbangkan goodness-of-fit dan kompleksitas.
Signifikansi Model - Variabel jumlah jam belajar signifikan meningkatkan peluang lulus, tetapi education numeric signifikan menurunkan peluang lulus. - Variabel gender juga signifikan pada taraf 5%.
Kesimpulan Praktis - Jumlah jam belajar lebih tinggi merupakan prediktor kuat untuk peluang lulus. - Gender juga berpotensi berpengaruh terhadap peluang lulus. - Model cukup baik dalam memprediksi dibandingkan model null.
Perbandingan Model
list(
AIC_Nominal = AIC(model_nominal),
AIC_Numeric = AIC(model_numeric)
)
## $AIC_Nominal
## [1] 679.3899
##
## $AIC_Numeric
## [1] 760.3701
Interpretasi: Semakin kecil nilai AIC, maka model tersebut semakin baik. Artinya, model nominal yang memperlakukan variabel ordinal menjadi nominal (dummy) lebih baik dibanding model numeric yang memperlakukan variabel ordinal menjadi rasio.
Goodness-of-Fit Model
nullmod <- glm(lulus ~ 1, data = data_ujian, family = binomial)
r2_nominal <- 1- (logLik(model_nominal)/logLik(nullmod))
r2_numeric <- 1- (logLik(model_numeric)/logLik(nullmod))
list(
McFadden_R2_Nominal = r2_nominal,
McFadden_R2_Numeric = r2_numeric
)
## $McFadden_R2_Nominal
## 'log Lik.' 0.1635136 (df=6)
##
## $McFadden_R2_Numeric
## 'log Lik.' 0.0570019 (df=4)
Interpretasi: McFadden’s R² mengukur goodness-of-fit. Semakin besar nilainya, semakin baik model memprediksi dibandingkan model null. Artinya, model nominal lebih baik dalam memprediksi peluang lulus dibandingkan model null.
Visualisasi Prediksi
sim_data_nominal <- data_ujian_nominal %>% mutate(predicted = predict(model_nominal, type = "response"))
sim_data_numeric <- data_ujian_numeric %>% mutate(predicted = predict(model_numeric, type = "response"))
# Plot untuk model nominal
library(ggplot2)
sim_data_nominal %>%
ggplot(aes(x = jumlah_jam_belajar, y = predicted, color = tingkat_pendidikan)) +
geom_point(alpha = 0.6) +
labs(
title = "Prediksi Probabilitas (Ordinal sebagai Nominal)",
x = "Jumlah Jam Belajar per Minggu",
y = "Probabilitas Lulus"
) +
theme_minimal()
sim_data_numeric %>%
ggplot(aes(x = jumlah_jam_belajar, y = predicted, color = as.factor(education_numeric))) +
geom_point(alpha = 0.6) +
labs(
title = "Prediksi Probabilitas (Ordinal sebagai Numeric)",
x = "Jumlah Jam Belajar per Minggu",
y = "Probabilitas Lulus"
) +
theme_minimal()
Interpretasi: Visualisasi hubungan antara jumlah jam belajar dan peluang lulus berdasarkan tingkat pendidikan, dengan dua pendekatan perlakuan ordinal.
Dalam membangun model regresi logistik, terdapat dua pendekatan utama yang sering digunakan, yaitu pendekatan Confirmatory dan Exploratory. Pendekatan ini memiliki perbedaan tujuan dan prosedur yang penting untuk dipahami oleh peneliti.
Pendekatan confirmatory menitikberatkan pada pengujian hipotesis yang sudah ditentukan sebelumnya. Model dibangun berdasarkan teori atau kajian literatur yang kuat, di mana variabel-variabel yang dimasukkan ke dalam model sudah dipilih berdasarkan hipotesis tertentu. Tujuan utamanya adalah untuk mengkonfirmasi apakah hubungan antar variabel sesuai dengan teori yang ada.
Kelebihan dari pendekatan ini adalah memberikan hasil yang lebih valid secara teoritis, serta meminimalkan risiko overfitting karena model tidak dibangun berdasarkan data semata. Namun, kekurangannya adalah kurang fleksibel jika terdapat variabel penting yang tidak dimasukkan dari awal.
Contoh
Berdasarkan kajian literatur yang sudah diketahui sebelumnya, didapatkan bahwa faktor usia, status merokok, dan indeks massa tubuh memengaruhi kesehatan jantung. Akan dibuat model regresi logistik berdasarkan variabel yang sudah diketahui tersebut untuk mengetahui apakah teori tersebut benar.
library(knitr)
library(dplyr)
library(ggplot2)
library(MASS)
library(DescTools)
set.seed(150)
jumlah_sampel <- 250
umur <- rnorm(jumlah_sampel)
status_merokok <- rbinom(jumlah_sampel, size = 1, prob = 0.5)
imtb <- rnorm(jumlah_sampel)
linear_pred <- -0.3 + 0.9 * umur - 1.1 * status_merokok + 0.5 * imtb
probabilitas <- 1 / (1 + exp(-linear_pred))
hasil <- rbinom(jumlah_sampel, size = 1, prob = probabilitas)
data_sintesis <- data.frame(
hasil = as.factor(hasil),
umur = umur,
merokok = status_merokok,
imtb = imtb
)
head(data_sintesis)
## hasil umur merokok imtb
## 1 0 -1.63230970 1 1.2240672
## 2 1 -0.06299626 0 0.1514091
## 3 1 -0.70544686 1 -0.2247016
## 4 0 -0.31417818 1 -1.0034459
## 5 0 -0.26694627 0 -0.7626163
## 6 1 0.15315947 0 -0.5497452
# Pemilihan Model
model_logistik <- glm(hasil ~ umur + merokok + imtb,
data = data_sintesis,
family = binomial())
summary(model_logistik)
##
## Call:
## glm(formula = hasil ~ umur + merokok + imtb, family = binomial(),
## data = data_sintesis)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1640 0.1957 -0.838 0.40205
## umur 0.8570 0.1725 4.967 6.79e-07 ***
## merokok -0.9352 0.2974 -3.144 0.00167 **
## imtb 0.7085 0.1711 4.141 3.45e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 328.94 on 249 degrees of freedom
## Residual deviance: 276.44 on 246 degrees of freedom
## AIC: 284.44
##
## Number of Fisher Scoring iterations: 4
Interpretasi: Seluruh variabel tersebut signifikan, artinya faktor-faktor tersebut benar-benar memengaruhi kesehatan jantung.
Pendekatan exploratory lebih fokus pada eksplorasi data untuk menemukan variabel-variabel yang signifikan secara empiris. Pendekatan ini biasanya digunakan ketika teori belum kuat atau ketika peneliti ingin menemukan pola baru dari data. Teknik-teknik seperti stepwise regression (baik forward maupun backward) sering digunakan untuk memilih variabel.
Pendekatan ini lebih fleksibel dan dapat mengungkap variabel prediktor yang sebelumnya tidak diperkirakan, tetapi berisiko menghasilkan model yang terlalu kompleks dan kurang generalisasi.
Metode stepwise regression adalah teknik pemilihan variabel dalam model regresi yang bertujuan menemukan kombinasi variabel terbaik untuk model prediksi. Ada tiga jenis utama:
Forward Selection: Mulai dari model tanpa variabel, kemudian variabel ditambahkan satu per satu berdasarkan kriteria signifikan sampai tidak ada variabel yang layak ditambahkan.
Backward Elimination: Mulai dari model penuh (semua variabel), kemudian variabel yang kurang signifikan dihapus satu per satu sampai hanya variabel signifikan tersisa.
Stepwise Selection (kedua arah): Kombinasi antara forward dan backward, dimana pada setiap langkah variabel dapat ditambahkan atau dihapus.
Metode ini membantu dalam membangun model yang sederhana namun kuat, menghindari overfitting dan meningkatkan interpretabilitas.
Contoh
Menggunakan contoh sebelumnya, kita ingin mengeksplorasi model untuk menemukan kombinasi variabel terbaik untuk model prediksi.
# Model full regresi logistik
model_logistik <- glm(hasil ~ umur + merokok + imtb,
data = data_sintesis,
family = binomial())
# Model null regresi logistik
null_model <- glm(hasil ~ 1, data = data_sintesis, family = binomial)
# Stepwise regression
step_forward <- step(null_model, direction = "forward", scope = formula(model_logistik), trace = FALSE)
step_backward <- step(model_logistik, direction = "backward", trace = FALSE)
step_both <- step(null_model, direction = "both", scope = formula(model_logistik), trace = FALSE)
AIC(model_logistik, step_forward, step_backward, step_both)
## df AIC
## model_logistik 4 284.4447
## step_forward 4 284.4447
## step_backward 4 284.4447
## step_both 4 284.4447
Interpretasi:
Semua model memiliki jumlah variabel yang sama, yaitu 4.
Nilai AIC pada semua model juga identik, yakni 284.4447.
Hal ini menunjukkan bahwa baik model penuh (model_logistik) maupun model hasil seleksi variabel dengan metode stepwise (forward, backward, dan kedua arah) memiliki performa yang sama menurut kriteria AIC.
ROC (Receiver Operating Characteristic) curve adalah grafik yang menampilkan performa model klasifikasi pada berbagai threshold probabilitas, dengan memplot true positive rate (sensitivitas) terhadap false positive rate (1 - spesifisitas). AUC (Area Under the Curve) mengukur keseluruhan performa model, dimana nilai AUC 1 menunjukkan model sempurna dan 0.5 menunjukkan model setara tebakan acak. Evaluasi ROC dan AUC sangat penting dalam regresi logistik untuk menilai kemampuan model membedakan kelas.
library(pROC)
## Warning: package 'pROC' was built under R version 4.3.3
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
pred_prob <- predict(step_both, type = "response")
roc_obj <- roc(data_sintesis$hasil, pred_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj, main = "Kurva ROC", col = "magenta")
auc(roc_obj)
## Area under the curve: 0.7616
Interpretasi: Berdasarkan perhitungan nilai AUC di atas, didapatkan nilai AUC sebesar 0.7616. Artinya, model memiliki kemampuan yang cukup baik dalam memisahkan dua kelas (negative dan positive).
Dalam regresi logistik, tidak ada R-squared seperti pada regresi linier. Oleh karena itu digunakan beberapa ukuran pseudo R-squared, seperti:
PseudoR2(step_both, which = c("CoxSnell", "Nagelkerke", "McFadden"))
## CoxSnell Nagelkerke McFadden
## 0.1894047 0.2588447 0.1595926
Interpretasi:
Cox & Snell R² (0.1894) Mengindikasikan bahwa model menjelaskan sekitar 18.94% variasi data.
Nagelkerke R² (0.2588) Nilai 0,2588 ini mengindikasikan bahwa model menjelaskan sekitar 25.88% variasi data. Artinya, model memiliki kemampuan penjelasan sedang terhadap variabilitas hasil.
McFadden R² (0.1596) Mengindikasikan bahwa model menjelaskan sekitar 15.96% variasi data.
# Mengubah probabilitas prediksi menjadi kelas biner
pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
# Membuat confusion matrix sederhana dengan fungsi table()
conf_matrix_alt <- table(Predicted = pred_class, Actual = as.numeric(as.character(data_sintesis$hasil)))
# Tampilkan confusion matrix
print(conf_matrix_alt)
## Actual
## Predicted 0 1
## 0 137 49
## 1 21 43
# Ekstrak nilai dari confusion matrix
TP <- conf_matrix_alt["1","1"]
TN <- conf_matrix_alt["0","0"]
FP <- conf_matrix_alt["1","0"]
FN <- conf_matrix_alt["0","1"]
# Hitung metrik
accuracy <- (TP + TN) / sum(conf_matrix_alt)
sensitivity <- TP / (TP + FN) # Recall
specificity <- TN / (TN + FP)
# Tampilkan metrik
list(
Accuracy = accuracy,
Sensitivity = sensitivity,
Specificity = specificity
)
## $Accuracy
## [1] 0.72
##
## $Sensitivity
## [1] 0.4673913
##
## $Specificity
## [1] 0.8670886
Interpretasi:
Model memiliki akurasi yang cukup baik yaitu sebesar 72%. Artinya, Model memprediksi dengan benar 72% dari keseluruhan data.
Nilai sensitivitas yang relatif rendah (46.7%) menunjukkan model kurang sensitif dalam mendeteksi kasus positif, sehingga ada risiko banyak kasus positif yang terlewat (false negatives).
Nilai spesifisitas yang tinggi (86.7%) menandakan model cukup baik dalam mengenali kasus negatif dan menghindari false positives.
Kesimpulan: Secara keseluruhan, model tersebut cukup baik karena memiliki akurasi sebesar 72% dan nilai AUC sebesar 76%. Namun, nilai sensitivitas sangat kecil padahal model tersebut diharapkan dapat lebih akurat dalam mendeteksi kasus positive (kesehatan jantung). Oleh karena itu, lebih baik menambahkan beberapa variabel dan mengevaluasi lagi model tersebut agar lebih akurat.
Bagian ini menyajikan cara membandingkan model regresi logistik menggunakan ukuran Deviance, AIC (Akaike Information Criterion), dan Likelihood-Ratio serta menjelaskan prinsip Parsimony dalam pemilihan model.
AIC
AIC adalah kriteria untuk memilih model terbaik dengan mempertimbangkan keseimbangan antara kompleksitas model dan goodness-of-fit.
Rumus AIC:
\[ \text{AIC} = -2 \times \log(L) + 2k \]
di mana:
AIC yang lebih rendah menunjukkan model yang lebih baik, karena memberikan penalti pada model yang terlalu kompleks.
Deviance
Deviance adalah ukuran seberapa baik model menjelaskan data dibandingkan dengan model sempurna (saturated model).
Rumus Deviance residual:
\[ D = -2 \times \log \left(\frac{L(\hat{\beta})}{L(\text{saturated})}\right) \]
Deviance yang lebih rendah menunjukkan model yang lebih baik fit-nya.
Contoh:
Seorang peneliti sedang melakukan studi untuk memprediksi apakah seseorang akan membeli produk tertentu (1 = membeli, 0 = tidak membeli) berdasarkan tiga variabel:
x1: Skor kepuasan pelanggan (numerik, berdistribusi normal), x2: Status keanggotaan pelanggan (biner, 1 = anggota, 0 = bukan anggota), x3: Skor aktivitas online pelanggan (numerik, berdistribusi normal). Peneliti mensimulasikan dataset dengan 300 observasi menggunakan model regresi logistik.
library(MASS)
library(broom)
## Warning: package 'broom' was built under R version 4.3.3
library(DescTools)
set.seed(123)
n <- 300
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
linear_predictor <- -0.5 + 0.9 * x1 + 1.1 * x2 - 0.7 * x3
prob <- 1 / (1 + exp(-linear_predictor))
y <- rbinom(n, 1, prob)
data <- data.frame(y = as.factor(y), x1, x2, x3)
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)
summary(model3)
##
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6365 0.1862 -3.418 0.000631 ***
## x1 0.7245 0.1486 4.874 1.09e-06 ***
## x2 1.3310 0.2649 5.025 5.03e-07 ***
## x3 -0.5607 0.1382 -4.056 4.98e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 415.77 on 299 degrees of freedom
## Residual deviance: 351.90 on 296 degrees of freedom
## AIC: 359.9
##
## Number of Fisher Scoring iterations: 3
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
## Model AIC Deviance
## 1 Model 1 400.2864 396.2864
## 2 Model 2 376.0173 370.0173
## 3 Model 3 359.9002 351.9002
Interpretasi: Semakin kecil nilai AIC dan Deviance, maka semakin baik juga model tersebut. Artinya, model 3 lebih baik daripada model 1 dan model 2 karena memiliki nilai AIC dan deviance terendah.
LRT digunakan untuk membandingkan dua model, yaitu model sederhana versus model yang lebih kompleks, guna menguji apakah penambahan parameter meningkatkan kualitas model secara signifikan.
Statistik uji LRT:
\[ \Lambda = -2 \left[ \log L(\text{model sederhana}) - \log L(\text{model kompleks}) \right] = D_{\text{sederhana}} - D_{\text{kompleks}} \]
Statistik \(\Lambda\) mengikuti distribusi \(\chi^2\) dengan derajat kebebasan sama dengan selisih jumlah parameter kedua model.
anova(model1, model2, test = "LRT")
## Analysis of Deviance Table
##
## Model 1: y ~ x1
## Model 2: y ~ x1 + x2
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 298 396.29
## 2 297 370.02 1 26.269 2.97e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Interpretasi: Nilai p-value = 0.0008909 (sangat kecil, jauh di bawah 0.05) menunjukkan bahwa penambahan variabel x2 secara signifikan meningkatkan kecocokan model terhadap data.
anova(model2, model3, test = "LRT")
## Analysis of Deviance Table
##
## Model 1: y ~ x1 + x2
## Model 2: y ~ x1 + x2 + x3
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 297 370.02
## 2 296 351.90 1 18.117 2.077e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Interpretasi: Nilai p-value = 1.153e-09 (sangat kecil, jauh di bawah 0.05) menunjukkan bahwa penambahan variabel x3 secara signifikan meningkatkan kecocokan model terhadap data.
Dalam pemodelan statistik, termasuk regresi logistik, prinsip parsimony (kadang disebut juga prinsip kesederhanaan atau Occam’s Razor) sangat penting untuk menghasilkan model yang efektif dan dapat diandalkan. Prinsip ini menekankan bahwa model yang paling sederhana yang mampu menjelaskan data dengan baik adalah model yang paling diinginkan.
Prinsip parsimony menyatakan bahwa:
“Antara dua model yang memiliki kemampuan prediksi yang hampir sama, model yang lebih sederhana lebih disukai.”
Hal ini berarti dalam konteks regresi logistik, tidak selalu menambahkan semua variabel prediktor ke dalam model adalah hal terbaik. Model yang terlalu kompleks dapat menyebabkan overfitting, yaitu model sangat cocok dengan data pelatihan tetapi buruk dalam memprediksi data baru.
Implikasi Prinsip Parsimony:
Menghindari overfitting Model dengan terlalu banyak variabel cenderung menangkap noise data, sehingga performa pada data baru menurun.
Interpretasi lebih mudah Model yang sederhana lebih mudah dipahami dan dijelaskan kepada pemangku kepentingan.
Efisiensi komputasi Model sederhana umumnya lebih cepat untuk diestimasi dan digunakan.
Model yang parsimonious dievaluasi dengan mempertimbangkan:
11.9.1 Sensitivitas dan Spesifitas
Sensitivitas: Kemampuan model mendeteksi kelas positif secara benar (True Positive Rate)
\[Sensitivity = \frac{TP}{TP + FN}\]
Spesifitas: Kemampuan model mendeteksi kelas negatif secara benar (True Negative Rate)
\[Specificity = \frac{TN}{TN + FP}\]
pred_prob <- predict(model3, type = "response")
# Mengubah probabilitas prediksi menjadi kelas biner
pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
# Membuat confusion matrix sederhana dengan fungsi table()
conf_matrix_alt <- table(Predicted = pred_class, Actual = as.numeric(as.character(data$y)))
# Tampilkan confusion matrix
print(conf_matrix_alt)
## Actual
## Predicted 0 1
## 0 92 46
## 1 55 107
# Ekstrak nilai dari confusion matrix
TP <- conf_matrix_alt["1","1"]
TN <- conf_matrix_alt["0","0"]
FP <- conf_matrix_alt["1","0"]
FN <- conf_matrix_alt["0","1"]
# Hitung metrik
accuracy <- (TP + TN) / sum(conf_matrix_alt)
sensitivity <- TP / (TP + FN) # Recall
specificity <- TN / (TN + FP)
# Tampilkan metrik
list(
Accuracy = accuracy,
Sensitivity = sensitivity,
Specificity = specificity
)
## $Accuracy
## [1] 0.6633333
##
## $Sensitivity
## [1] 0.6993464
##
## $Specificity
## [1] 0.6258503
Interpretasi: - Model memiliki akurasi yang cukup baik yaitu sebesar 70%. Artinya, Model memprediksi dengan benar 70% dari keseluruhan data. - Nilai sensitivitas yang relatif rendah (72%) menunjukkan model cukup sensitif dalam mendeteksi kasus positif. - Nilai spesifisitas yang tinggi (69%) menandakan model kurang baik dalam mengenali kasus negatif dan menghindari false positives.
Kurva ROC Kurva ROC merupakan alat visual untuk menilai performa model klasifikasi biner. Kurva ini menggambarkan hubungan antara True Positive Rate (Sensitivity) dan False Positive Rate (1 - Specificity) pada berbagai tingkat threshold klasifikasi. 1. Definisi
Sumbu Y: Sensitivity = True Positive Rate = \(\frac{TP}{TP + FN}\) Sumbu X: 1 - Specificity = False Positive Rate = \(\frac{FP}{FP + TN}\) Garis diagonal (dari pojok kiri bawah ke kanan atas) mencerminkan performa acak (random guess). Kurva yang mendekati sudut kiri atas menandakan performa model klasifikasi yang lebih baik.
Sensitivity meningkat Specificity menurun
Sebaliknya, ketika cut-off dinaikkan, model menjadi lebih konservatif:
Sensitivity menurun Specificity meningkat
Naik secara vertikal dengan cepat hingga sensitivity mencapai 1 BergerContained secara horizontal menuju 1 - specificity = 1 Area under the curve (AUC) mendekati nilai 1
AUC = 0.5: model tidak lebih baik dari tebak acak AUC > 0.7: model cukup baik AUC > 0.9: model sangat baik AUC dikenal juga sebagai concordance index, yaitu probabilitas bahwa model memberikan nilai skor probabilitas yang lebih tinggi untuk kasus positif daripada kasus negatif.
Untuk membandingkan performa beberapa model klasifikasi Untuk memilih threshold (cut-off) optimal berdasarkan kebutuhan aplikasi (misalnya: lebih penting menghindari false negative atau false positive?)
Kurva ROC dapat dibuat menggunakan package pROC:
library(pROC)
pred_prob <- predict(model3, type = "response")
roc_obj <- roc(data$y, pred_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj, main = "Kurva ROC", col = "purple")
auc(roc_obj)
## Area under the curve: 0.7545
Interpretasi: Berdasarkan perhitungan nilai AUC di atas, didapatkan nilai AUC sebesar 0.7832. Artinya, model memiliki kemampuan yang cukup baik dalam memisahkan dua kelas (negative dan positive).
Untuk memilih threshold terbaik, kita bisa mengevaluasi sensitivitas dan spesifisitas pada berbagai cut-off.
thresholds <- seq(0.1, 0.9, by = 0.05)
results <- data.frame(Threshold = thresholds)
results$Sensitivity <- sapply(thresholds, function(t) {
pred_class <- ifelse(pred_prob >= t, 1, 0)
cm <- table(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_prob >= 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.10 1.00000000 0.02040816
## 2 0.15 0.99346405 0.08163265
## 3 0.20 0.98039216 0.17006803
## 4 0.25 0.94117647 0.23129252
## 5 0.30 0.91503268 0.36734694
## 6 0.35 0.86928105 0.48299320
## 7 0.40 0.80392157 0.55102041
## 8 0.45 0.76470588 0.58503401
## 9 0.50 0.69934641 0.62585034
## 10 0.55 0.66013072 0.70068027
## 11 0.60 0.54248366 0.78231293
## 12 0.65 0.48366013 0.86394558
## 13 0.70 0.35947712 0.91156463
## 14 0.75 0.28104575 0.95238095
## 15 0.80 0.16993464 0.98639456
## 16 0.85 0.08496732 0.99319728
## 17 0.90 0.02614379 0.99319728
Interpretasi:
Ketika threshold menurun (misalnya dari 0.55 ke 0.10), Sensitivity meningkat (dari 0.6624 ke 1.0000), tetapi Specificity menurun (dari 0.7273 ke 0.0699). Ini menunjukkan trade-off klasik: model menjadi lebih sensitif dalam mendeteksi kasus positif (kurang melewatkan kasus positif), tetapi kurang spesifik (lebih banyak salah mengklasifikasikan kasus negatif sebagai positif).
Sebaliknya, ketika threshold naik, Sensitivity menurun dan Specificity meningkat, menunjukkan model menjadi lebih konservatif dan lebih baik dalam mengidentifikasi kasus negatif.
Karena peneliti ingin memprediksi apakah seseorang akan membeli produknya atau tidak, maka peneliti lebih baik meminimalkan false negative (FN). False negative artinya model memprediksi seseorang tidak akan membeli (0), tetapi sebenarnya mereka membeli (1). Ini berarti peluang penjualan hilang karena pelanggan yang tertarik tidak ditargetkan dengan kampanye pemasaran atau penawaran. Dengan meminimalkan FN, perusahaan dapat memastikan lebih banyak pelanggan potensial ditargetkan, sehingga hal tersebut dapat menangkap sebanyak mungkin peluang penjualan.
Precision-Recall Curve (PR Curve) merupakan instrumen untuk menilai efektivitas model klasifikasi, terutama ketika menghadapi dataset dengan distribusi kelas yang tidak merata (class imbalance). Alat ini lebih informatif dibandingkan Kurva ROC dalam situasi di mana kelas minoritas memiliki signifikansi tinggi, karena PR Curve fokus pada performa terhadap kelas positif.
Presisi: Persentase prediksi positif yang ternyata benar-benar positif, mengukur akurasi dari prediksi positif.
\[Precision = \frac{TP}{TP + FP}\]
Recall: Persentase kasus positif yang berhasil diidentifikasi dengan benar oleh model, mencerminkan kemampuan model untuk menangkap semua kasus positif.
\[Recall = \frac{TP}{TP + FN}\]
Kurva PR lebih disukai daripada Kurva ROC ketika dataset memiliki ketidakseimbangan kelas yang parah, karena ROC Curve dapat memberikan gambaran yang terlalu optimis dengan menekankan true negative rate yang kurang relevan dalam kasus tersebut.
| Aspek | ROC Curve | Precision-Recall Curve |
|---|---|---|
| Fokus | Semua kelas | Kelas positif saja |
| Kuat di | Data seimbang | Data tidak seimbang |
| Sumbu Y | Sensitivitas (Recall) | Precision |
| Sumbu X | 1 - Spesifisitas | Recall |
Luas kurva (AUPRC) mendekati 1 berarti model sangat baik. Baseline AUPRC = prevalensi kelas positif dalam data.
Kurva PR sering digunakan dalam aplikasi seperti deteksi penyakit langka atau klasifikasi teks (misalnya, mendeteksi spam), di mana kelas positif (penyakit atau spam) jauh lebih sedikit dibandingkan kelas negatif, tetapi memiliki dampak besar jika terlewat.
library(PRROC)
## Warning: package 'PRROC' was built under R version 4.3.3
## Loading required package: rlang
set.seed(456) # Ubah seed untuk menghasilkan data berbeda
x1 <- rnorm(200)
x2 <- rbinom(200, 1, 0.4) # Ubah probabilitas x2 menjadi 0.4
x3 <- rnorm(200)
lin_pred <- -2 + 1.2 * x1 - 0.5 * x2 + 0.8 * x3 # Ubah koefisien: -2, 1.2, -0.5, 0.8
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(200, 1, p)
data <- data.frame(y = y, x1, x2, x3)
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
prob <- predict(model, type = "response")
pr <- pr.curve(scores.class0 = prob[data$y == 1],
scores.class1 = prob[data$y == 0],
curve = TRUE)
plot(pr)
Bagian ini akan menjelaskan perhitungan R-squared pada regresi logistik yaitu \(R^2_{\text{CoxandSnell}}\) dan \(R^2_{\text{McFadden}}\).
Definisi: \(R^2_{\text{CoxandSnell}}\) adalah ukuran pseudo R-squared yang mengukur proporsi variasi dalam data yang dijelaskan oleh model logistik, berdasarkan pendekatan likelihood ratio. Ini menyesuaikan dengan fakta bahwa model logistik tidak memiliki R-squared tradisional.
Rumus:
\(R^2_{\text{Cox and Snell}} = 1 - \left( \frac{L_0}{L_M} \right)^{2/n}\)
Dengan:
\(L_0\): Likelihood model null (tanpa prediktor)
\(L_M\): Likelihood model penuh
\(n\): Jumlah observasi
Interpretasi: Nilai berkisar antara 0 dan kurang dari 1 (biasanya maksimum sekitar 0.75-0.8 pada dataset besar). Nilai yang lebih tinggi menunjukkan model yang lebih baik dalam menjelaskan data.
Definisi: \(R^2_{\text{McFadden}}\) adalah pseudo R-squared lain yang mengukur perbaikan model logistik dibandingkan dengan model null berdasakan rasio log-likelihood.
Rumus:
\(R^2_{\text{McFadden}} = 1 - \frac{\log(L_M)}{\log(L_0)}\)
Dengan:
\(L_0\): Likelihood model null (tanpa prediktor)
\(L_M\): Likelihood model penuh
Interpretasi: Nilai berkisar antara 0 dan 1. Nilai kecil (misalnya, < 0.2) menunjukkan model yang moderat, sedangkan nilai > 0.4 dianggap sangat baik, meskipun ini tergantung pada konteks dataset.
\(R^2_{\text{CoxandSnell}}\) cenderung memberikan nilai yang lebih konservatif dibandingkan \(R^2_{\text{McFadden}}\) dan sering digunakan untuk membandingkan model dengan jumlah parameter yang berbeda.
\(R^2_{\text{McFadden}}\) lebih sederhana dan sering dipilih untuk menilai kebaikan model secara relatif.
set.seed(123)
n <- 300
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
linear_predictor <- -0.5 + 0.9 * x1 + 1.1 * x2 - 0.7 * x3
prob <- 1 / (1 + exp(-linear_predictor))
y <- rbinom(n, 1, prob)
data <- data.frame(y = as.factor(y), x1, x2, x3)
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
model_null <- glm(y ~ 1, data = data, family = binomial)
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
## R2_Cox_Snell R2_McFadden
## 1 0.191758 0.1536147
Perhitungan Otomatis dengan Package Tambahan Menggunakan pscl
library(pscl)
## Warning: package 'pscl' was built under R version 4.3.3
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
pR2(model)
## fitting null model for pseudo-r2
## llh llhNull G2 McFadden r2ML r2CU
## -175.9500943 -207.8841502 63.8681118 0.1536147 0.1917580 0.2557114
library(DescTools)
PseudoR2(model, which = "all")
## McFadden McFaddenAdj CoxSnell Nagelkerke AldrichNelson
## 0.1536147 0.1343732 0.1917580 0.2557114 0.1755254
## VeallZimmermann Efron McKelveyZavoina Tjur AIC
## 0.3021768 0.1943572 0.2569516 0.1957186 359.9001886
## BIC logLik logLik0 G2
## 374.7153185 -175.9500943 -207.8841502 63.8681118
Interpretasi:
Cox & Snell R² (0.192) Mengindikasikan bahwa model menjelaskan sekitar 19,2% variasi data.
McFadden R² (0.154) Mengindikasikan bahwa model menjelaskan sekitar 15.4% variasi data.
Berdasarkan nilai R-squared, dapat dikatakan bahwa model tersebut memiliki kekuatan yang kurang baik dalam prediktif.
Distribusi multinomial merupakan pengembangan dari distribusi binomial yang digunakan untuk menangani lebih dari dua kategori.
Apabila \(X_1, X_2, \ldots, X_k\) mewakili jumlah kemunculan pada masing-masing dari \(k\) kategori, maka probabilitasnya dapat ditulis 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 syarat \(\sum_{i=1}^k x_i = n\) dan \(\sum_{i=1}^k p_i = 1\).
Contoh Kasus Distribusi Multinomial
Seorang peneliti melakukan survei terhadap 100 orang untuk mengetahui preferensi mereka terhadap tiga merek minuman: A, B, dan C. Berdasarkan data sebelumnya, diketahui bahwa 40% memilih merek A, 35% memilih merek B, dan 25% memilih merek C. Dari 100 orang yang disurvei, asumsikan hasilnya mengikuti distribusi multinomial. Jika ternyata 45 orang memilih merek A, 30 orang memilih merek B, dan 25 orang memilih merek C, hitung probabilitas kombinasi ini terjadi.
Jawab:
Jumlah total percobaan: \(n = 100\)
Jumlah kejadian untuk setiap kategori: \(X_1 = 45\) (merek A) \(X_2 = 30\) (merek B) \(X_3 = 25\) (merek C)
Probabilitas untuk setiap kategori: \(p_1 = 0.40\) (merek A) \(p_2 = 0.35\) (merek B) \(p_3 = 0.25\) (merek C)
Verifikasi: \(\sum_{i=1}^{3} p_i = 0.40 + 0.35 + 0.25 = 1\), yang sesuai.
Probabilitas dari distribusi multinomial diberikan oleh: \[P(X_1 = x_1, X_2 = x_2, X_3 = x_3) = \frac{n!}{x_1! x_2! x_3!} p_1^{x_1} p_2^{x_2} p_3^{x_3}\] dengan \(\sum_{i=1}^{3} x_i = n\) dan \(\sum_{i=1}^{3} p_i = 1\).
Substitusikan nilai ke dalam rumus:
\(n! = 100!\) \(x_1! = 45!\), \(x_2! = 30!\), \(x_3! = 25!\) \(p_1^{x_1} = (0.40)^{45}\), \(p_2^{x_2} = (0.35)^{30}\), \(p_3^{x_3} = (0.25)^{25}\)
\[P(X_1 = 45, X_2 = 30, X_3 = 25) = \frac{100!}{45! \cdot 30! \cdot 25!} \cdot (0.40)^{45} \cdot (0.35)^{30} \cdot (0.25)^{25}\]
probs <- c(0.40, 0.35, 0.25)
counts <- c(45, 30, 25)
probability <- dmultinom(counts, prob = probs, log = FALSE)
cat("Probabilitas:", probability, "\n")
## Probabilitas: 0.004376503
Hasilnya (dihitung dengan R): Probabilitas ≈ 0,0044.
Kesimpulan: probabilitas bahwa 45 orang memilih merek A, 30 orang memilih merek B, dan 25 orang memilih merek C dari 100 orang adalah sekitar 0,0044, yang menunjukkan kombinasi ini cukup jarang terjadi berdasarkan probabilitas yang diberikan, tetapi masih mungkin dalam distribusi multinomial.
Model ini diterapkan untuk menganalisis hubungan antara variabel respons dengan lebih dari dua kategori dan satu atau lebih variabel prediktor. Jika variabel respons \(Y\) memiliki \(K\) kategori, dan kategori ke-\(K\) dipilih sebagai kategori acuan (kategori dasar), maka model logistik untuk kategori \(j\) dapat dirumuskan 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\).
###Baseline-category logit model
Baseline-category logit model merupakan pendekatan regresi logistik untuk variabel respons dengan lebih dari dua kategori (kategori nominal). Dalam model ini, satu kategori ditetapkan sebagai kategori acuan (kategori dasar), dan probabilitas kategori lainnya dibandingkan dengan kategori acuan dalam bentuk logit:
\[\log\left(\frac{\pi_j}{\pi_c}\right), \quad j = 1, \ldots, c-1\]
dengan:
\(\pi_j\): Probabilitas respons berada pada kategori \(j\) \(\pi_c\): Probabilitas respons berada pada kategori acuan (kategori dasar)
Dengan demikian, model ini menghasilkan \((c-1)\) fungsi logit.
Model regresi logistik dengan satu prediktor \(x\) memiliki bentuk umum:
Misalkan \(Y \in \{1, 2, 3\}\) dengan kategori 3 sebagai baseline:
Kategori 1 vs. baseline:
\[ \log\left(\frac{\pi_1}{\pi_3}\right) = \alpha_1 + \beta_1 x \]
Kategori 2 vs. baseline:
\[ \log\left(\frac{\pi_2}{\pi_3}\right) = \alpha_2 + \beta_2 x \]
Perbedaan log-odds antara kategori 1 dan 2:
\[ \log\left(\frac{\pi_1}{\pi_2}\right) = \log\left(\frac{\pi_1 / \pi_3}{\pi_2 / \pi_3}\right) = (\alpha_1 + \beta_1 x) - (\alpha_2 + \beta_2 x) = (\alpha_1 - \alpha_2) + (\beta_1 - \beta_2) x \]
Model Baseline-Category Logit:
Estimasi parameter dilakukan dengan metode maximum likelihood menggunakan algoritma iteratif seperti Newton-Raphson. Fungsi log-likelihood didefinisikan sebagai:
\[ \ell(\beta) = \sum_{i=1}^{n} \sum_{j=1}^{K} y_{ij} \log(\pi_{ij}) \]
di mana \(\pi_{ij} = P(Y_i = j | x_i)\) adalah probabilitas respons \(Y_i\) pada kategori \(j\) untuk observasi \(i\), dan \(y_{ij} = 1\) jika \(Y_i = j\), serta \(y_{ij} = 0\) sebaliknya.
Algoritma Newton-Raphson memperbarui parameter berdasarkan gradien dan Hessian dari fungsi log-likelihood.
Sebuah perusahaan teknologi ingin memahami faktor-faktor yang memengaruhi preferensi pengguna terhadap jenis layanan cloud computing, yaitu Shared Hosting, VPS, atau Dedicated Server.
Perusahaan mensurvei 200 pengguna dan mengumpulkan data berikut:
Tujuan: Mengetahui pengaruh kelompok usia, bidang pekerjaan, dan pengalaman terhadap pilihan layanan cloud computing.
Simulasi Data
set.seed(125) # Mengubah seed untuk variasi
n <- 200 # Mengubah jumlah pengguna dari 200 menjadi 220
Job <- sample(c("Developer", "Data Analyst", "Manager"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 38, sd = 8)) # Mengubah mean dari 35 menjadi 38 dan sd dari 7 menjadi 8
Experience <- round(pmax(rnorm(n, mean = 6, sd = 2.5), 0)) # Mengubah mean dari 7 menjadi 6 dan sd dari 3 menjadi 2.5
# Simulasikan Service berdasarkan probabilitas berbeda per Job
Service <- sapply(Job, function(job) {
if (job == "Developer") {
sample(c("Shared Hosting", "VPS", "Dedicated Server"), size = 1, prob = c(0.5, 0.3, 0.2))
} else if (job == "Data Analyst") {
sample(c("Shared Hosting", "VPS", "Dedicated Server"), size = 1, prob = c(0.3, 0.4, 0.3))
} else {
sample(c("Shared Hosting", "VPS", "Dedicated Server"), size = 1, prob = c(0.2, 0.5, 0.3))
}
})
df <- data.frame(Service = factor(Service), Age, Job = factor(Job), Experience)
df$Service <- relevel(df$Service, ref = "Shared Hosting") # baseline
head(df)
## Service Age Job Experience
## 1 Shared Hosting 51 Data Analyst 7
## 2 VPS 46 Data Analyst 12
## 3 VPS 47 Manager 3
## 4 Dedicated Server 33 Manager 4
## 5 Shared Hosting 42 Developer 4
## 6 Shared Hosting 47 Developer 3
Estimasi Model
library(nnet)
## Warning: package 'nnet' was built under R version 4.3.3
model_mnlogit <- multinom(Service ~ Age + Job + Experience, data = df)
## # weights: 18 (10 variable)
## initial value 219.722458
## iter 10 value 192.322223
## final value 191.613655
## converged
Interpretasi:
summary(model_mnlogit)
## Call:
## multinom(formula = Service ~ Age + Job + Experience, data = df)
##
## Coefficients:
## (Intercept) Age JobDeveloper JobManager Experience
## Dedicated Server 2.4163440 -0.05955408 -1.541560 -0.01482368 0.002464637
## VPS 0.8036673 -0.03658591 -2.018879 0.67465109 0.128165800
##
## Std. Errors:
## (Intercept) Age JobDeveloper JobManager Experience
## Dedicated Server 1.076159 0.02431405 0.4739776 0.4606144 0.08025350
## VPS 1.057146 0.02338692 0.5175491 0.4284842 0.07851647
##
## Residual Deviance: 383.2273
## AIC: 403.2273
Interpretasi:
Koefisien menunjukkan perubahan log-odds relatif terhadap baseline: Dedicated Server:
(Intercept): 2.4163440Age: -0.05955408 → Peningkatan 1 tahun usia mengurangi
log-odds sebesar 0.0596.JobDeveloper: -1.541560 → Developer memiliki log-odds
1.5416 lebih rendah.JobManager: -0.01482368 → Manager memiliki log-odds
0.0148 lebih rendah.Experience: 0.002464637 → Peningkatan 1 tahun
pengalaman menambah log-odds 0.0025. VPS:(Intercept): 0.8036673Age: -0.03658591 → Peningkatan 1 tahun usia mengurangi
log-odds sebesar 0.0366.JobDeveloper: -2.018879 → Developer memiliki log-odds
2.0189 lebih rendah.JobManager: 0.67465109 → Manager memiliki log-odds
0.6747 lebih tinggi.Experience: 0.128165800 → Peningkatan 1 tahun
pengalaman menambah log-odds 0.1282.Standar error mengukur ketidakpastian estimasi. Nilai kecil
menunjukkan estimasi presisi, seperti Experience
(0.08025350 untuk Dedicated Server, 0.07851647 untuk VPS).
Residual Deviance: 383.2273 → Menunjukkan penyimpangan model.
AIC: 403.2273 → Digunakan untuk membandingkan model; nilai lebih rendah lebih baik.
Kesimpulan:
Nilai P-Value dan Interpretasi
z <- summary(model_mnlogit)$coefficients / summary(model_mnlogit)$standard.errors
pval <- 2 * (1- pnorm(abs(z)))
round(pval, 4)
## (Intercept) Age JobDeveloper JobManager Experience
## Dedicated Server 0.0247 0.0143 0.0011 0.9743 0.9755
## VPS 0.4471 0.1177 0.0001 0.1154 0.1026
Interpretasi:
Nilai p-value kecil (<0.05) menunjukkan variabel tersebut signifikan memengaruhi preferensi pilihan layanan cloud computing. Berdasarkan perhitungan di atas, dapat disimpulkan bahwa:
Dedicated Server: Usia dan JobDeveloper berpengaruh nyata, sedangkan JobManager dan Experience tidak.
VPS: Hanya JobDeveloper berpengaruh nyata, sedangkan Age, JobManager, dan Experience tidak signifikan.
Prediksi dan Validasi
df$Predicted <- predict(model_mnlogit)
table(Predicted = df$Predicted, Actual = df$Service)
## Actual
## Predicted Shared Hosting Dedicated Server VPS
## Shared Hosting 50 14 13
## Dedicated Server 7 16 11
## VPS 22 23 44
Interpretasi:
Kesimpulan:
Model tersebut kurang akurat dalam memprediksi preferensi pilihan layanan cloud computing. Selain itu, banyak variabel yang tidak signifikan dalam memengaruhi preferensi pilihan layanan cloud computing. Model tersebut dapat diperbaiki dengan menambah variabel atau mengganti variabel yang signifikan.
Regresi logistik ordinal digunakan ketika variabel dependen \(Y\) memiliki sifat ordinal, yaitu kategori yang memiliki urutan tertentu, namun jarak antar kategori tidak harus sama. Contohnya adalah tingkat kepuasan pelanggan yang dapat dikelompokkan menjadi: Rendah, Sedang, dan Tinggi.
Model ini berbeda dengan:
Model yang sering digunakan dalam regresi logistik ordinal adalah Cumulative Logit Model, yang mengasumsikan proportional odds. Model ini memodelkan logaritma odds kumulatif sebagai fungsi linear dari variabel prediktor:
\[ \log \left( \frac{P(Y \leq j)}{P(Y > j)} \right) = \alpha_j + \beta x \]
Di mana:
Untuk \(c\) kategori, model ini membangun \(c-1\) persamaan logit kumulatif.
Koefisien \(\beta\) dalam model regresi logistik multinomial menunjukkan efek prediktor \(x\) terhadap probabilitas memilih kategori tertentu dibandingkan kategori lain, baik yang lebih rendah maupun lebih tinggi:
Odds ratio mengukur rasio odds antar kategori untuk setiap peningkatan satu unit pada prediktor \(x\), dengan rumus:
\[ \text{OR} = e^{\beta} \]
Model regresi logistik berurutan yang sering digunakan adalah pendekatan Cumulative Logit Model dengan asumsi Proportional Odds, juga disebut hipotesis garis sejajar. Hipotesis paralelisme menyatakan bahwa koefisien regresi (\(\beta\)) untuk prediktor konsisten di seluruh kategori kumulatif variabel respons.
Model ditulis sebagai:
\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j + \beta x \]
untuk \(j = 1, \ldots, c-1\), di mana:
Karakteristik Model: - Intersep (\(\alpha_j\)) berbeda antar batas kategori. - Koefisien \(\beta\) seragam di seluruh fungsi logit kumulatif. - Visualisasi logit kumulatif menunjukkan kemiringan paralel, hanya berbeda posisi intersep.
Regresi Logistik Oerdinal harus memenuhi asumsi paralelisme. Bebereapa risiko jika asumsi paralelisme tidak terpenuhi adalah sebagai berikut: - Efek prediktor bervariasi per batas kategori. - Model Cumulative Logit tidak valid. - Alternatif: Adjacent-category logit dan Continuation-ratio (sequential) logit
Untuk memeriksa validitas asumsi, dapat digunakan: - Likelihood Ratio Test antara model proportional dan non-proportional
Uji ini membandingkan model dengan asumsi proportional odds (terbatas) dan model tanpa asumsi tersebut (tidak terbatas).
Metode: Hitung log-likelihood dari kedua model, lalu gunakan statistik \(-2 \times (\text{log-likelihood terbatas} - \text{log-likelihood tidak terbatas})\) yang mengikuti distribusi chi-square. Interpretasi: Jika p-value < 0.05, asumsi paralelisme ditolak
Uji ini, dari package brant di R, memeriksa konsistensi koefisien antar batas kategori.
Metode: Bandingkan model ordinal dengan model di mana koefisien bervariasi per batas, menggunakan statistik chi-square. Interpretasi: Jika p-value prediktor < 0.05, asumsi paralelisme tidak terpenuhi
Contoh
Sebuah perusahaan pengembang aplikasi ingin menganalisis data untuk mengevaluasi tingkat kepuasan pengguna terhadap layanan aplikasi mereka (1: Tidak Puas, 2: Cukup Puas, 3: Sangat Puas) berdasarkan kecepatan respons dan frekuensi penggunaan. Perusahaan mengumpulkan data dari 250 pengguna.
set.seed(456)
n <- 250
response_time <- round(runif(n, 1, 10))
usage_frequency <- cut(runif(n), breaks = c(-Inf, 0.4, 0.7, Inf), labels = c("Rendah", "Sedang", "Tinggi"))
satisfaction <- cut(runif(n), breaks = c(-Inf, 0.3, 0.6, Inf), labels = c("Tidak Puas", "Cukup Puas", "Sangat Puas"))
df <- data.frame(satisfaction = factor(satisfaction, ordered = TRUE), response_time, usage_frequency = factor(usage_frequency))
head(df)
## satisfaction response_time usage_frequency
## 1 Cukup Puas 2 Rendah
## 2 Tidak Puas 3 Rendah
## 3 Cukup Puas 8 Sedang
## 4 Sangat Puas 9 Sedang
## 5 Cukup Puas 8 Rendah
## 6 Sangat Puas 4 Rendah
Uji Paralelisme dengan Brant Test
Hipotesis Nol (\(H_0\)): Koefisien regresi (\(\beta\)) untuk prediktor sama di seluruh batas kategori, sehingga asumsi paralelisme terpenuhi.
Hipotesis Alternatif (\(H_1\)): Koefisien regresi untuk prediktor berbeda antar batas kategori, sehingga asumsi paralelisme tidak terpenuhi.
alpha: 0,05
library(MASS)
library(brant)
## Warning: package 'brant' was built under R version 4.3.3
model_ord <- polr(satisfaction ~ response_time + usage_frequency, data = df, Hess = TRUE)
brant(model_ord)
## ----------------------------------------------------
## Test for X2 df probability
## ----------------------------------------------------
## Omnibus 1.43 3 0.7
## response_time 0.65 1 0.42
## usage_frequencySedang 0.17 1 0.68
## usage_frequencyTinggi 0.23 1 0.63
## ----------------------------------------------------
##
## H0: Parallel Regression Assumption holds
Tolak H0 jika p-value < 0,05
Keputusan
Omnibus:
p-value 0.69 (> 0.05) menunjukkan tidak ada bukti pelanggaran asumsi
paralelisme secara keseluruhan.
response time:
p-value 0.24 (> 0.05), terima H0.
usage frequency sedang:
p-value 0.91 (> 0.05), terima H0.
usage frequency tinggi:
p-value 0.79 (> 0.05), terima H0.
Kesimpulan
Tidak ada pelanggaran asumsi paralelisme untuk semua prediktor dan secara keseluruhan (p-value > 0.05).
Model Cumulative Logit Model valid untuk analisis lebih lanjut.
Estimasi Model Ordinal
library(MASS)
model_ord <- polr(satisfaction ~ response_time + usage_frequency, data = df, Hess = TRUE)
summary(model_ord)
## Call:
## polr(formula = satisfaction ~ response_time + usage_frequency,
## data = df, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## response_time 0.008925 0.04554 0.19600
## usage_frequencySedang -0.025716 0.29368 -0.08756
## usage_frequencyTinggi -0.041766 0.27968 -0.14933
##
## Intercepts:
## Value Std. Error t value
## Tidak Puas|Cukup Puas -1.0561 0.3407 -3.0997
## Cukup Puas|Sangat Puas 0.1603 0.3336 0.4806
##
## Residual Deviance: 529.4895
## AIC: 539.4895
(ctable <- coef(summary(model_ord)))
## Value Std. Error t value
## response_time 0.008925439 0.04553887 0.19599608
## usage_frequencySedang -0.025715696 0.29368353 -0.08756261
## usage_frequencyTinggi -0.041765558 0.27968412 -0.14933118
## Tidak Puas|Cukup Puas -1.056124513 0.34071490 -3.09973095
## Cukup Puas|Sangat Puas 0.160295188 0.33356028 0.48055837
Interpretasi Koefisien Koefisien Prediktor:
response_time: 0.06120
Setiap peningkatan 1 detik pada response time meningkatkan log-odds
kepuasan yang lebih tinggi sebesar 0.06120.
Odds Ratio: \(( e^{0.06120} \approx 1.063
)\), odds naik 6.3%.
usage_frequencySedang: 0.07264
Frekuensi Sedang meningkatkan log-odds kepuasan sebesar 0.07264
dibandingkan Rendah.
Odds Ratio: \(( e^{0.07264} \approx 1.075
)\), odds naik 7.5%.
usage_frequencyTinggi: 0.28645
Frekuensi Tinggi meningkatkan log-odds sebesar 0.28645 dibandingkan
Rendah.
Odds Ratio: \(( e^{0.28645} \approx 1.332
)\), odds naik 33.2%.
Intersep (Thresholds):
Tidak Puas|Cukup Puas: -0.5789 (SE 0.3208, t = -1.8045)
Log-odds dasar untuk transisi dari Tidak Puas ke Cukup Puas atau lebih
tinggi adalah -0.5789 saat response time = 0 dan usage frequency =
Rendah. Nilai negatif berarti odds untuk Cukup Puas atau lebih tinggi
turun sebesar 0.5789 dari nol, menunjukkan kecenderungan awal untuk
tetap di Tidak Puas.
Odds: \(( e^{-0.5789} \approx 0.560
)\), odds 44% lebih kecil.
Cukup Puas|Sangat Puas: 0.7997 (SE 0.3227, t = 2.4779)
Log-odds dasar untuk transisi dari Cukup Puas ke Sangat Puas adalah
0.7997, menunjukkan odds untuk Sangat Puas naik sebesar 0.7997 dari nol,
mencerminkan kecenderungan signifikan untuk beralih ke kategori lebih
tinggi.
Odds: \(( e^{0.7997} \approx 2.225 )\),
odds naik 122.5%.
Residual Deviance dan AIC: - Residual Deviance: 537.3842 - AIC:
547.3842
- Digunakan untuk membandingkan model alternatif.
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = round(p, 4)))
## Value Std. Error t value p value
## response_time 0.008925439 0.04553887 0.19599608 0.8446
## usage_frequencySedang -0.025715696 0.29368353 -0.08756261 0.9302
## usage_frequencyTinggi -0.041765558 0.27968412 -0.14933118 0.8813
## Tidak Puas|Cukup Puas -1.056124513 0.34071490 -3.09973095 0.0019
## Cukup Puas|Sangat Puas 0.160295188 0.33356028 0.48055837 0.6308
Interpretasi: Koefisien Prediktor: - response time: t = 1.4059224, p
= 0.1597
Tidak signifikan (p > 0.05), waktu respons tidak memengaruhi kepuasan
secara nyata. - usage frequency sedang: t = 0.2487967, p = 0.8035
Tidak signifikan (p > 0.05), frekuensi Sedang tidak berpengaruh. -
usage frequency tinggi: t = 1.0183161, p = 0.3085
Tidak signifikan (p > 0.05), frekuensi Tinggi juga tidak
berpengaruh.
Intersep (Thresholds): - Tidak Puas|Cukup Puas: t = -1.8044774, p =
0.0712
Tidak signifikan (p > 0.05), transisi ini kurang jelas. - Cukup
Puas|Sangat Puas: t = 2.4778556, p = 0.0132
Signifikan (p < 0.05), transisi ini menunjukkan perbedaan nyata.
Kesimpulan: Karena tidak ada prediktor yang signifikan, sebaiknya ditambahkan lagi prediktor lain atau mengganti prediktor pada model tersebut.
newdata <- data.frame(response_time = 5:9, usage_frequency = "Sedang")
predict(model_ord, newdata = newdata, type = "probs")
## Tidak Puas Cukup Puas Sangat Puas
## 1 0.2544466 0.2808406 0.4647128
## 2 0.2527571 0.2803091 0.4669338
## 3 0.2510750 0.2797689 0.4691560
## 4 0.2494005 0.2792200 0.4713795
## 5 0.2477334 0.2786625 0.4736041
Interpretasi: Probabilitas Sangat Puas naik dari 39.63% (5 detik) ke 45.60% (9 detik), sementara Tidak Puas dan Cukup Puas turun, menunjukkan waktu respons lebih lama dikaitkan dengan kepuasan lebih tinggi pada kategori tertinggi.
Log Linear Model adalah pendekatan statistik untuk menganalisis hubungan antar variabel kategorik dalam tabel kontingensi, terutama fokus pada frekuensi atau jumlah pengamatan. Model ini termasuk dalam kelas Generalized Linear Model (GLM) dan memodelkan log frekuensi yang diharapkan (\(\log \mu_{ij\ldots}\)) berdasarkan efek utama variabel dan interaksi di antaranya. Berbeda dengan regresi logistik yang berfokus pada variabel respons biner atau ordinal, Log Linear Model lebih cocok untuk mengeksplorasi struktur asosiasi antar variabel tanpa variabel dependen yang jelas, sering digunakan dalam analisis data kategorik multidimensi.
Model log-linier mengasumsikan log frekuensi yang diharapkan dalam sel tabel kontingensi dinyatakan sebagai kombinasi linier efek utama dan interaksi. Untuk tabel dua arah dengan variabel \(A\) dan \(B\), model ditulis sebagai:
\[ \log \mu_{ij} = \lambda + \lambda_i^A + \lambda_j^B + \lambda_{ij}^{AB} \]
di mana:
- \(\mu_{ij}\) adalah frekuensi yang diharapkan untuk kombinasi kategori \(i\) dari \(A\) dan \(j\) dari \(B\),
- \(\lambda\) adalah efek konstanta,
- \(\lambda_i^A\) dan \(\lambda_j^B\) adalah efek utama variabel \(A\) dan \(B\),
- \(\lambda_{ij}^{AB}\) adalah efek interaksi antara \(A\) dan \(B\).
Jika tidak ada interaksi (\(\lambda_{ij}^{AB} = 0\)), model mengasumsikan kemandirian antara \(A\) dan \(B\), yang dapat diuji dengan statistik chi-square atau likelihood ratio test.
Log Linear Model sering digunakan di sosiologi, kedokteran, dan ilmu sosial untuk menganalisis data seperti hubungan antara jenis kelamin, status perkawinan, dan tingkat pendapatan, atau antara faktor risiko dan penyakit. Misalnya, dalam epidemiologi, model ini dapat mengevaluasi interaksi antara merokok dan paparan asbes terhadap risiko kanker paru. Model ini juga fleksibel untuk tabel dengan lebih dari dua dimensi, memungkinkan pemodelan interaksi tingkat tinggi.
Parameter diestimasi dengan metode maksimum likelihood. Pemilihan model dilakukan secara bertahap, dari model sederhana (hanya efek utama) hingga model dengan interaksi. Uji likelihood ratio (\(G^2\)) atau chi-square (\(X^2\)) digunakan untuk membandingkan model, di mana:
\[ G^2 = 2 \sum (O_{ij} \log \frac{O_{ij}}{E_{ij}}) \]
dengan \(O_{ij}\) sebagai frekuensi teramati dan \(E_{ij}\) sebagai frekuensi yang diharapkan. Nilai p menentukan apakah model sederhana cukup atau perlu interaksi tambahan.
Dalam analisis data kategorik, terdapat beberapa pendekatan diantaranya yaitu, tabel kontingensi, Log Linear Model, dan Model Regresi Logistik. Ketiga pendekatan tersebut memiliki perbedaan masing-masing yaitu:
| Aspek | Tabel Kontingensi | Model Loglinier | Regresi Logistik |
|---|---|---|---|
| Tujuan | Deskripsi frekuensi | Deteksi asosiasi (simetris) | Prediksi probabilitas |
| Variabel dependen | Tidak ada | Tidak ada (simetris) | Ada (eksplisit) |
| Distribusi | Tidak diasumsikan | Poisson (frekuensi sel) | Binomial (probabilitas) |
| Bentuk Model | Tidak ada | GLM: \(\log(\mu) \sim efek\) | GLM: \(logit(p) \sim prediktor\) |
| Cocok untuk | Eksplorasi awal | Tabel > 2 variabel | Studi prediktif |
| Asumsi Model | Tidak memerlukan asumsi distribusi atau model statistik | Mengasumsikan distribusi Poisson, dapat diperluas dengan pendekatan Bayesian | Mengasumsikan distribusi binomial, dengan linieritas dalam logit |
| Fleksibilitas Dimensi | Terbatas pada data dua atau tiga arah tanpa pemodelan formal | Fleksibel untuk tabel multidimensi, mendukung interaksi tingkat tinggi | Fokus pada hubungan prediktor-respon dengan dimensi sederhana |
| Output | Frekuensi atau persentase mentah | Estimasi efek dan pengujian asosiasi | Probabilitas atau odds ratio untuk prediksi |
Dalam analisis log-linear model, terdapat dua jenis model yang sering digunakan untuk menganalisis hubungan antar variabel kategorikal dalam data kontingensi: model saturasi (saturated model) dan model independen (independent model). Kedua model ini memiliki pendekatan yang berbeda dalam memodelkan interaksi antar variabel, dengan model saturasi mempertimbangkan semua interaksi yang mungkin, sementara model independen mengasumsikan bahwa variabel-variabel tersebut tidak saling berinteraksi. Pemilihan model yang tepat bergantung pada struktur data dan tujuan analisis.
Model-model ini sering dibandingkan menggunakan uji likelihood ratio untuk menentukan model mana yang paling sesuai dengan data yang dianalisis.
Model saturasi adalah model log-linear yang paling kompleks, yang mencakup semua efek utama dan interaksi antar variabel dalam tabel kontingensi. Dalam model ini, diasumsikan bahwa semua hubungan yang mungkin antara kategori dalam variabel-variabel yang terlibat ada dan harus diperhitungkan. Dengan kata lain, model saturasi menganggap bahwa tidak ada asumsi independensi antar variabel dan bahwa setiap kombinasi interaksi antara kategori dari setiap variabel adalah signifikan.
Model saturasi untuk tiga variabel kategorikal dapat dirumuskan sebagai berikut:
\[ \log(E_{ijk}) = \mu + \lambda_i + \lambda_j + \lambda_k + \lambda_{ij} + \lambda_{ik} + \lambda_{jk} + \lambda_{ijk} \]
Dimana:
\(E_{ijk}\) adalah ekspektasi frekuensi untuk kategori \(i\), \(j\), dan \(k\).
\(\mu\) adalah rata-rata atau intercept model.
\(\lambda_i, \lambda_j, \lambda_k\) adalah efek utama dari masing-masing kategori pada variabel \(i\), \(j\), dan \(k\).
\(\lambda_{ij}, \lambda_{ik}, \lambda_{jk}\) adalah efek interaksi antara dua variabel.
\(\lambda_{ijk}\) adalah efek interaksi tiga variabel.
Model ini memungkinkan peneliti untuk menangkap semua interaksi yang mungkin antara kategori dalam data, tetapi bisa menjadi sangat rumit ketika jumlah kategori atau variabel bertambah, memerlukan banyak parameter untuk diestimasi, dan sering kali membutuhkan ukuran sampel yang besar agar hasilnya stabil.
Kelebihan Model Saturasi: - Kemampuan menangkap semua hubungan antar variabel: Model ini mempertimbangkan semua interaksi yang mungkin, memberikan gambaran yang sangat rinci tentang hubungan antar kategori. - Dapat menangani data yang lebih kompleks: Sangat berguna ketika ada alasan kuat untuk percaya bahwa interaksi antar variabel mempengaruhi hasil yang diamati.
Kekurangan Model Saturasi: - Model yang rumit: Semakin banyak variabel yang dianalisis, semakin banyak parameter yang harus diestimasi, membuat model menjadi lebih rumit dan sulit diinterpretasikan. - Memerlukan data yang besar: Karena model ini memiliki banyak parameter, dibutuhkan ukuran sampel yang besar untuk mendapatkan estimasi yang akurat dan stabil.
Model independen, di sisi lain, adalah model yang lebih sederhana, yang mengasumsikan bahwa tidak ada interaksi antara variabel-variabel yang dianalisis. Dalam model ini, hanya efek utama dari masing-masing variabel yang diperhitungkan, tanpa mempertimbangkan interaksi antara kategori variabel. Model ini mengasumsikan bahwa setiap variabel bersifat independen satu sama lain.
Model independen untuk tiga variabel kategorikal dapat dituliskan sebagai berikut:
\[ \log(E_{ijk}) = \mu + \lambda_i + \lambda_j + \lambda_k \]
Dimana:
\(E_{ijk}\) adalah ekspektasi frekuensi untuk kategori \(i\), \(j\), dan \(k\).
\(\mu\) adalah rata-rata atau intercept model.
\(\lambda_i, \lambda_j, \lambda_k\) adalah efek utama dari masing-masing kategori pada variabel \(i\), \(j\), dan \(k\).
Model ini lebih sederhana dan lebih mudah diestimasi karena hanya melibatkan efek utama dari setiap variabel. Namun, model ini mengabaikan interaksi yang mungkin ada antara variabel-variabel, yang dapat menyebabkan hilangnya informasi penting dalam beberapa kasus.
Kelebihan Model Independen: - Kesederhanaan: Model ini mudah dipahami dan diestimasi, serta memiliki jumlah parameter yang lebih sedikit. - Cepat dan efisien: Karena hanya melibatkan efek utama, estimasi parameter bisa dilakukan lebih cepat dan memerlukan ukuran sampel yang lebih kecil.
Kekurangan Model Independen: - Mengabaikan interaksi: Model ini tidak memperhitungkan interaksi antar variabel, yang bisa mengarah pada kesimpulan yang tidak akurat jika interaksi antar variabel memang signifikan. - Kurang fleksibel: Model ini hanya cocok untuk data yang tidak menunjukkan interaksi signifikan antar variabel. Jika ada interaksi yang kuat, model ini akan kurang efektif.
Perbandingan Model Saturasi dan Model Independen:
| Aspek | Model Saturasi | Model Independen |
|---|---|---|
| Kompleksitas | Sangat kompleks, mempertimbangkan semua interaksi | Sederhana, hanya mempertimbangkan efek utama |
| Jumlah Parameter | Banyak parameter (semua interaksi) | Hanya mempertimbangkan efek utama |
| Asumsi | Tidak ada asumsi independensi antar variabel | Asumsi independensi antar variabel |
| Kelebihan | Dapat menangkap semua hubungan antar variabel | Mudah diestimasi dan diinterpretasi |
| Kekurangan | Memerlukan banyak data dan parameter, sulit diinterpretasi | Bisa mengabaikan interaksi penting antar variabel |
Untuk memilih antara model saturasi dan model independen, kita dapat menggunakan uji likelihood ratio (LRT). Uji ini membandingkan kedua model berdasarkan seberapa baik masing-masing model menjelaskan data yang diamati.
Dalam uji likelihood ratio, kita menguji dua hipotesis sebagai berikut:
Contoh
Sebuah survei dilakukan untuk menganalisis preferensi pekerja terhadap metode kerja jarak jauh (remote work) berdasarkan status pekerjaan mereka. Variabel yang diamati adalah Status Pekerjaan (Penuh Waktu atau Paruh Waktu) dan Preferensi Remote (Setuju atau Tidak Setuju). Data frekuensi dari survei disajikan dalam tabel kontingensi 2x2 berikut:
| Status Pekerjaan | Preferensi Setuju | Preferensi Tidak Setuju |
|---|---|---|
| Penuh Waktu | 45 | 30 |
| Paruh Waktu | 20 | 55 |
Lakukan analisis dengan saturated model dan independent model, kemudian uji apakah ada interaksi yang signifikan antara Status Pekerjaan dan Preferensi Remote menggunakan uji likelihood ratio test.
Jawab:
Tabel Kontingensi dan Model Loglinier
# Membuat tabel kontingensi 2x2
tabel <- matrix(c(45, 30, 20, 55), nrow = 2,
dimnames = list(Status_Pekerjaan = c("Penuh Waktu", "Paruh Waktu"),
Preferensi_Remote = c("Setuju", "Tidak Setuju")))
print(tabel)
## Preferensi_Remote
## Status_Pekerjaan Setuju Tidak Setuju
## Penuh Waktu 45 20
## Paruh Waktu 30 55
Rumus Model: Model log-linier untuk tabel ini dapat ditulis sebagai:
\[ \log(\mu_{ij}) = \mu + \lambda_i^P + \lambda_j^R + \lambda_{ij}^{PR} \]
di mana \(\mu\) adalah efek rata-rata, \(\lambda_i^P\) adalah efek utama status pekerjaan, \(\lambda_j^R\) adalah efek utama preferensi remote, dan \(\lambda_{ij}^{PR}\) adalah efek interaksi antara status pekerjaan dan preferensi remote.
Model Saturated
library(MASS)
model_saturated <- loglm(~ Status_Pekerjaan * Preferensi_Remote, data = tabel)
summary(model_saturated)
## Formula:
## ~Status_Pekerjaan * Preferensi_Remote
## attr(,"variables")
## list(Status_Pekerjaan, Preferensi_Remote)
## attr(,"factors")
## Status_Pekerjaan Preferensi_Remote
## Status_Pekerjaan 1 0
## Preferensi_Remote 0 1
## Status_Pekerjaan:Preferensi_Remote
## Status_Pekerjaan 1
## Preferensi_Remote 1
## attr(,"term.labels")
## [1] "Status_Pekerjaan" "Preferensi_Remote"
## [3] "Status_Pekerjaan:Preferensi_Remote"
## 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
Interpretasi Elemen Spesifik 1. Angka
1 0 1 0 1 1 1 2 dan Artinya - factors: - Nilai:
1 0 1 (efek utama dan interaksi untuk
Status_Pekerjaan), 0 1 1 (efek utama dan
interaksi untuk Preferensi_Remote). - Makna: Mendefinisikan
struktur model dengan efek utama dan interaksi. - order [1 1 2]: -
Urutan efek: 1 (efek utama Status_Pekerjaan,
Preferensi_Remote), 2 (interaksi). - Makna:
Mencerminkan kompleksitas istilah, dengan interaksi memiliki order lebih
tinggi.
Model Independent
model_indep <- loglm(~ Status_Pekerjaan + Preferensi_Remote, data = tabel)
summary(model_indep)
## Formula:
## ~Status_Pekerjaan + Preferensi_Remote
## attr(,"variables")
## list(Status_Pekerjaan, Preferensi_Remote)
## attr(,"factors")
## Status_Pekerjaan Preferensi_Remote
## Status_Pekerjaan 1 0
## Preferensi_Remote 0 1
## attr(,"term.labels")
## [1] "Status_Pekerjaan" "Preferensi_Remote"
## 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 17.33050 1 3.141041e-05
## Pearson 16.96833 1 3.800865e-05
Interpretasi Output 1. Formula dan Struktur Model: -
~ Status_Pekerjaan + Preferensi_Remote: Hanya efek utama,
tanpa interaksi. - attr(,"term.labels") dan
attr(,"order") [1 1]: Dua efek utama (order
1). - attr(,"intercept") [1]: Adanya intersep
sebagai baseline. - attr(,"response") [0]:
Tidak ada variabel respons eksplisit.
Odds Ratio
Odds Ratio dihitung untuk mengukur asosiasi antara Status Pekerjaan dan Preferensi Remote menggunakan rumus:
\[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \]
di mana:
\(n_{11} = 45\) (Penuh Waktu, Setuju),
\(n_{12} = 30\) (Penuh Waktu, Tidak Setuju),
\(n_{21} = 20\) (Paruh Waktu, Setuju),
\(n_{22} = 55\) (Paruh Waktu, Tidak Setuju).
# Menghitung Odds Ratio secara manual
n11 <- tabel[1, 1] # Penuh Waktu, Setuju
n12 <- tabel[1, 2] # Penuh Waktu, Tidak Setuju
n21 <- tabel[2, 1] # Paruh Waktu, Setuju
n22 <- tabel[2, 2] # Paruh Waktu, Tidak Setuju
odds_ratio <- (n11 * n22) / (n12 * n21)
cat("Odds Ratio (manual):", odds_ratio, "\n")
## Odds Ratio (manual): 4.125
Interpretasi \(OR > 1\) menunjukkan asosiasi positif. Peluang Penuh Waktu setuju dengan remote work 4.125 kali lebih besar dibandingkan Paruh Waktu.
Estimasi Parameter
logOR <- log((tabel[1,1] * tabel[2,2]) / (tabel[1,2] * tabel[2,1]))
logOR
## [1] 1.417066
Log OR dapat digunakan sebagai estimasi awal parameter interaksi dalam model jenuh.
Likelihood Ratio Test
anova(model_indep, model_saturated)
## LR tests for hierarchical log-linear models
##
## Model 1:
## ~Status_Pekerjaan + Preferensi_Remote
## Model 2:
## ~Status_Pekerjaan * Preferensi_Remote
##
## Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1 17.3305 1
## Model 2 0.0000 0 17.3305 1 3e-05
## Saturated 0.0000 0 0.0000 0 1e+00
Interpretasi Signifikansi Interaksi: \(P = 3e-05\) menolak hipotesis nol kemandirian, menunjukkan adanya interaksi signifikan. Artinya, model saturated lebih baik.
Model log-linear yang melibatkan tiga variabel kategorik, seperti X, Y, dan Z, memberikan pendekatan untuk menganalisis interaksi antara variabel-variabel tersebut. Model log-linear dapat dikembangkan dalam berbagai bentuk tergantung pada sejauh mana interaksi yang ingin dimasukkan dalam model. Model-model ini digunakan untuk memahami hubungan antara variabel yang saling terkait dalam tabel kontingensi, dengan tujuan utama untuk mengeksplorasi efek utama dan interaksi antara variabel.
Berikut adalah beberapa model log-linear yang sering digunakan dalam analisis tabel tiga arah:
Model saturated adalah model yang mencakup seluruh interaksi yang mungkin terjadi antara variabel-variabel yang ada. Dalam model ini, tidak ada batasan mengenai jenis interaksi yang dapat terjadi, sehingga mencakup interaksi antara dua variabel serta interaksi tiga arah antara ketiga variabel.
Model ini dapat dituliskan sebagai berikut:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} + \lambda_{ijk}^{XYZ} \]
Dengan model ini, kita dapat menganalisis secara komprehensif bagaimana setiap variabel dan kombinasi interaksinya mempengaruhi hasil yang diamati.
Model homogen menyederhanakan model saturated dengan hanya mempertimbangkan interaksi dua arah antara variabel-variabel yang ada, tanpa memperhitungkan interaksi tiga arah. Ini memberikan model yang lebih sederhana dan dapat berguna dalam beberapa kasus di mana interaksi tiga arah dianggap tidak signifikan atau tidak relevan.
Model ini dituliskan sebagai berikut:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \]
Dalam hal ini, hanya interaksi dua variabel yang dimasukkan, yang lebih sederhana dibandingkan model saturated.
Model conditional berfokus pada hubungan antar variabel dengan mempertimbangkan satu variabel sebagai kondisi untuk variabel lainnya. Model ini menguji efek interaksi antar variabel dengan mengontrol variabel tertentu.
Conditional pada X:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} \]
Model ini memuat interaksi antara X dengan Y dan antara X dengan Z.
Conditional pada Y:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{jk}^{YZ} \]
Model ini memuat interaksi antara Y dengan X dan antara Y dengan Z.
Conditional pada Z:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \]
Model ini memuat interaksi antara Z dengan X dan antara Z dengan Y.
Model ini mengasumsikan adanya independensi antara pasangan variabel, dan hanya mempertimbangkan interaksi antara pasangan variabel yang terkait. Ada tiga bentuk independensi yang dapat diuji dalam model ini:
Independensi antara X & Y:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} \]
Independensi antara X & Z:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ik}^{XZ} \]
Independensi antara Y & Z:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{jk}^{YZ} \]
Model ini memberikan pandangan yang lebih terbatas tentang hubungan antara variabel dengan mengasumsikan bahwa dua variabel tertentu tidak berinteraksi satu sama lain.
Dalam model ini, tidak ada interaksi antar variabel yang dipertimbangkan. Model ini hanya memperhitungkan efek utama dari setiap variabel tanpa mempertimbangkan kombinasi interaksi apa pun di antara mereka.
Model ini dituliskan sebagai berikut:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z \]
Model ini paling sederhana dan sering digunakan ketika interaksi antar variabel dianggap tidak penting atau tidak relevan.
Pengujian Interaksi dalam Model Log-Linear Tiga Arah
Pengujian interaksi dalam analisis model log-linear tiga arah bertujuan untuk mengidentifikasi apakah terdapat hubungan interaksi yang signifikan antara variabel-variabel yang ada. Proses pengujian ini dilakukan secara bertahap, dimulai dari pengujian interaksi yang paling kompleks (tiga arah) hingga yang lebih sederhana (dua arah atau tanpa interaksi). Berikut adalah tahapan pengujian interaksi yang umum dilakukan pada model log-linear dengan tiga variabel kategorik (X, Y, Z):
Setiap tahap pengujian bertujuan untuk mengevaluasi kesesuaian model dan menentukan struktur interaksi mana yang paling tepat sesuai dengan data yang diamati. Pengujian interaksi ini membantu dalam menentukan model yang paling parsimonious (sederhana namun memadai) untuk menggambarkan hubungan antara variabel.
Contoh Soal
Tabel berikut menyajikan data dari survei yang dilakukan pada kelompok umur, jenis pekerjaan, dan sikap terhadap kebijakan peningkatan pajak untuk produk tembakau. Susun dan interpretasikan model log-linear paling sederhana (paling parsimonious) untuk data ini. Jelaskan proses yang Anda lakukan dalam menentukan model terbaik serta asosiasi apa saja yang teridentifikasi. Tunjukkan juga bagaimana nilai yang diprediksi dari model menggambarkan asosiasi tersebut.
| Kelompok Umur | Jenis Pekerjaan | Mendukung | Menolak | Total |
|---|---|---|---|---|
| 18-25 tahun | Mahasiswa | 80 | 40 | 120 |
| 18-25 tahun | Pekerja | 70 | 50 | 120 |
| 18-25 tahun | Total | 150 | 90 | 240 |
| 26-40 tahun | Mahasiswa | 60 | 30 | 90 |
| 26-40 tahun | Pekerja | 90 | 60 | 150 |
| 26-40 tahun | Total | 150 | 90 | 240 |
| 41-60 tahun | Mahasiswa | 50 | 20 | 70 |
| 41-60 tahun | Pekerja | 80 | 40 | 120 |
| 41-60 tahun | Total | 130 | 60 | 190 |
Keterangan:
Jawab:
library("epitools")
library("DescTools")
library("lawstat")
## Warning: package 'lawstat' was built under R version 4.3.3
Input data
# Input data sesuai tabel praktikum
z.age <- factor(rep(c("18-25", "26-40", "41-60"), each = 4))
x.job <- factor(rep(c("Student", "Worker"), each = 2, times = 3))
y.tax <- factor(rep(c("Support", "Oppose"), times = 6))
counts <- c(80, 40, 70, 50, 60, 30, 90, 60, 50, 20, 80, 40)
data <- data.frame(
Age_Group = z.age,
Job_Type = x.job,
Tax_Support = y.tax,
Frequency = counts
)
data
## Age_Group Job_Type Tax_Support Frequency
## 1 18-25 Student Support 80
## 2 18-25 Student Oppose 40
## 3 18-25 Worker Support 70
## 4 18-25 Worker Oppose 50
## 5 26-40 Student Support 60
## 6 26-40 Student Oppose 30
## 7 26-40 Worker Support 90
## 8 26-40 Worker Oppose 60
## 9 41-60 Student Support 50
## 10 41-60 Student Oppose 20
## 11 41-60 Worker Support 80
## 12 41-60 Worker Oppose 40
Membentuk tabel kontingensi 3 arah
# Membuat tabel kontingensi 3 arah
table3d <- xtabs(Frequency ~ Age_Group + Job_Type + Tax_Support, data = data)
ftable(table3d)
## Tax_Support Oppose Support
## Age_Group Job_Type
## 18-25 Student 40 80
## Worker 50 70
## 26-40 Student 30 60
## Worker 60 90
## 41-60 Student 20 50
## Worker 40 80
UJI MODEL INTERAKSI TIGA ARAH (SATURATED VS HOMOGENOUS)
# Penentuan kategori reference
x.job <- relevel(x.job, ref = "Worker")
y.tax <- relevel(y.tax, ref = "Oppose")
z.age <- relevel(z.age, ref = "26-40")
Model Saturated memasukkan semua interaksi hingga tiga arah:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} + \lambda_{ijk}^{XYZ} \]
# Model Saturated
model_saturated <- glm(Frequency ~ Job_Type + Tax_Support + Age_Group +
Job_Type*Tax_Support + Job_Type*Age_Group + Tax_Support*Age_Group +
Job_Type*Tax_Support*Age_Group,
data = data, family = poisson(link = "log"))
summary(model_saturated)
##
## Call:
## glm(formula = Frequency ~ Job_Type + Tax_Support + Age_Group +
## Job_Type * Tax_Support + Job_Type * Age_Group + Tax_Support *
## Age_Group + Job_Type * Tax_Support * Age_Group, family = poisson(link = "log"),
## data = data)
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 3.689e+00 1.581e-01 23.331
## Job_TypeWorker 2.231e-01 2.121e-01 1.052
## Tax_SupportSupport 6.931e-01 1.936e-01 3.579
## Age_Group26-40 -2.877e-01 2.415e-01 -1.191
## Age_Group41-60 -6.931e-01 2.739e-01 -2.531
## Job_TypeWorker:Tax_SupportSupport -3.567e-01 2.679e-01 -1.331
## Job_TypeWorker:Age_Group26-40 4.700e-01 3.082e-01 1.525
## Job_TypeWorker:Age_Group41-60 4.700e-01 3.464e-01 1.357
## Tax_SupportSupport:Age_Group26-40 1.532e-15 2.958e-01 0.000
## Tax_SupportSupport:Age_Group41-60 2.231e-01 3.279e-01 0.681
## Job_TypeWorker:Tax_SupportSupport:Age_Group26-40 6.899e-02 3.867e-01 0.178
## Job_TypeWorker:Tax_SupportSupport:Age_Group41-60 1.335e-01 4.234e-01 0.315
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## Job_TypeWorker 0.292841
## Tax_SupportSupport 0.000344 ***
## Age_Group26-40 0.233608
## Age_Group41-60 0.011373 *
## Job_TypeWorker:Tax_SupportSupport 0.183113
## Job_TypeWorker:Age_Group26-40 0.127286
## Job_TypeWorker:Age_Group41-60 0.174850
## Tax_SupportSupport:Age_Group26-40 1.000000
## Tax_SupportSupport:Age_Group41-60 0.496136
## Job_TypeWorker:Tax_SupportSupport:Age_Group26-40 0.858410
## Job_TypeWorker:Tax_SupportSupport:Age_Group41-60 0.752486
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 9.6196e+01 on 11 degrees of freedom
## Residual deviance: 2.6645e-14 on 0 degrees of freedom
## AIC: 93.397
##
## Number of Fisher Scoring iterations: 3
# Menghitung eksponensial koefisien untuk interpretasi
exp(model_saturated$coefficients)
## (Intercept)
## 40.000000
## Job_TypeWorker
## 1.250000
## Tax_SupportSupport
## 2.000000
## Age_Group26-40
## 0.750000
## Age_Group41-60
## 0.500000
## Job_TypeWorker:Tax_SupportSupport
## 0.700000
## Job_TypeWorker:Age_Group26-40
## 1.600000
## Job_TypeWorker:Age_Group41-60
## 1.600000
## Tax_SupportSupport:Age_Group26-40
## 1.000000
## Tax_SupportSupport:Age_Group41-60
## 1.250000
## Job_TypeWorker:Tax_SupportSupport:Age_Group26-40
## 1.071429
## Job_TypeWorker:Tax_SupportSupport:Age_Group41-60
## 1.142857
Interpretasi Koefisien
(Intercept): Rata-rata log jumlah kasus untuk
kategori referensi (Mahasiswa, Menolak Pajak, 18-25 tahun)
adalah 3.69
(atau \(\mu \approx 40\)).
Job_TypeWorker: Pekerja memiliki expected count sekitar 1.25 kali lebih tinggi dibandingkan Mahasiswa dalam kategori referensi, namun tidak signifikan (p = 0.2928).
Tax_SupportSupport: Mereka yang mendukung pajak tembakau memiliki expected count sekitar 2 kali lipat dibandingkan yang menolak, dengan signifikansi yang kuat (p = 0.000344).
Age_Group26-40: Kelompok umur 26-40 tahun memiliki expected count sekitar 0.75 kali lebih rendah dibandingkan dengan kelompok referensi 18-25 tahun, namun tidak signifikan (p = 0.2336).
Age_Group41-60: Kelompok umur 41-60 tahun memiliki expected count sekitar 0.5 kali lebih rendah dibandingkan dengan kelompok referensi 18-25 tahun, dan signifikan (p = 0.01137).
Interaksi dua arah:
Job_TypeWorker:Tax_SupportSupport: Pekerja yang mendukung pajak tembakau memiliki expected count sekitar 0.7 kali lebih rendah dibandingkan Mahasiswa yang mendukung pajak tembakau, namun tidak signifikan (p = 0.1831).
Job_TypeWorker:Age_Group26-40: Pekerja yang berada di kelompok umur 26-40 tahun memiliki expected count sekitar 1.6 kali lebih tinggi dibandingkan Mahasiswa dalam kelompok umur 26-40 tahun, namun tidak signifikan (p = 0.1273).
Job_TypeWorker:Age_Group41-60: Pekerja yang berada di kelompok umur 41-60 tahun memiliki expected count sekitar 1.6 kali lebih tinggi dibandingkan Mahasiswa dalam kelompok umur 41-60 tahun, namun tidak signifikan (p = 0.1749).
Interaksi tiga arah:
Job_TypeWorker:Tax_SupportSupport:Age_Group26-40: Pekerja yang mendukung pajak tembakau dalam kelompok umur 26-40 tahun memiliki expected count sekitar 1.07 kali lebih tinggi dibandingkan Mahasiswa yang mendukung pajak tembakau dalam kelompok umur 26-40 tahun, namun tidak signifikan (p = 0.8584).
Job_TypeWorker:Tax_SupportSupport:Age_Group41-60: Pekerja yang mendukung pajak tembakau dalam kelompok umur 41-60 tahun memiliki expected count sekitar 1.14 kali lebih tinggi dibandingkan Mahasiswa yang mendukung pajak tembakau dalam kelompok umur 41-60 tahun, namun tidak signifikan (p = 0.7525).
Goodness-of-Fit
Residual deviance ≈ 0 menunjukkan bahwa model saturated sangat fit terhadap data (seluruh variasi data dijelaskan oleh model).
AIC = 93.397 menunjukkan nilai yang baik untuk model ini, dan dapat digunakan untuk perbandingan dengan model-model lain yang lebih sederhana.
Kesimpulan
Model saturated ini sangat fit dengan data, dengan residual deviance yang sangat rendah menunjukkan bahwa model ini memadai dalam menggambarkan data.
Efek utama yang paling signifikan adalah:
Sikap mendukung pajak tembakau: expected count untuk yang mendukung pajak tembakau adalah 2 kali lebih tinggi dibandingkan yang menolak (signifikan, p = 0.000344).
Kelompok umur 41-60 tahun: expected count untuk kelompok umur 41-60 tahun adalah 0.5 kali lebih rendah dibandingkan dengan kelompok umur 18-25 tahun (signifikan, p = 0.01137).
Tidak ditemukan bukti kuat adanya interaksi dua atau tiga arah yang signifikan, dengan sebagian besar interaksi tidak signifikan (p > 0.05).
Model yang lebih sederhana (tanpa interaksi tiga arah) perlu dipertimbangkan untuk model final yang lebih parsimonious.
Catatan Interpretasi:
Nilai exp(coef) menunjukkan rasio
ekspektasi dibandingkan dengan kategori referensi.
Efek positif → menaikkan expected
count;
Efek negatif → menurunkan expected
count.
Koefisien signifikan memiliki nilai p-value < 0.05.
Model Homogenous
Model log-linear homogenous memasukkan semua efek utama dan semua interaksi dua arah, tanpa interaksi tiga arah. Model homogenus:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \]
Uji Hipotesis: Apakah Ada Interaksi Tiga Arah? (Saturated vs Homogenous)
Pengujian ini bertujuan untuk membandingkan model saturated dan model homogenous dengan menggunakan residual deviance dari kedua model tersebut.
Hipotesis
Hitung Selisih Deviance
# Model Homogen
model_homogenous <- glm(Frequency ~ Job_Type + Tax_Support + Age_Group +
Job_Type*Tax_Support + Job_Type*Age_Group + Tax_Support*Age_Group,
family = poisson(link = "log"), data = data)
summary(model_homogenous)
##
## Call:
## glm(formula = Frequency ~ Job_Type + Tax_Support + Age_Group +
## Job_Type * Tax_Support + Job_Type * Age_Group + Tax_Support *
## Age_Group, family = poisson(link = "log"), data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.70920 0.13916 26.654 < 2e-16 ***
## Job_TypeWorker 0.18627 0.16627 1.120 0.262603
## Tax_SupportSupport 0.66251 0.15919 4.162 3.16e-05 ***
## Age_Group26-40 -0.31235 0.18818 -1.660 0.096942 .
## Age_Group41-60 -0.74862 0.20911 -3.580 0.000344 ***
## Job_TypeWorker:Tax_SupportSupport -0.29782 0.16664 -1.787 0.073904 .
## Job_TypeWorker:Age_Group26-40 0.51340 0.18608 2.759 0.005796 **
## Job_TypeWorker:Age_Group41-60 0.55916 0.19908 2.809 0.004973 **
## Tax_SupportSupport:Age_Group26-40 0.03716 0.19017 0.195 0.845096
## Tax_SupportSupport:Age_Group41-60 0.30266 0.20706 1.462 0.143823
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 96.19552 on 11 degrees of freedom
## Residual deviance: 0.10131 on 2 degrees of freedom
## AIC: 89.498
##
## Number of Fisher Scoring iterations: 3
# Deviance antar model
Deviance.model <- model_homogenous$deviance - model_saturated$deviance
Deviance.model
## [1] 0.1013108
Hitung derajat bebas
# Derajat bebas (df) antar model
derajat.bebas <- model_homogenous$df.residual - model_saturated$df.residual
derajat.bebas
## [1] 2
Chi-Square Tabel (α = 0.05)
# Menghitung nilai chi-square dari tabel dengan derajat bebas
chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel
## [1] 5.991465
Keputusan uji
# Membandingkan selisih deviance dengan nilai chi-square untuk membuat keputusan
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0", "Tolak H0")
Keputusan
## [1] "Terima H0"
Interpretasi: Pada taraf nyata 5%, belum cukup bukti untuk menolak H0. Artinya tidak ada interaksi tiga arah antara jenis kelamin, fundamentalisme, dan pendapat mengenai hukuman mati.
Rangkuman pengujian
Hipotesis
Hitung Selisih Deviance
Dari hasil output model, kita menghitung selisih deviance antara model homogenous dan model saturated:
\[ \Delta \text{Deviance} = \text{Deviance model homogenous} - \text{Deviance model saturated} \]
Hasilnya adalah:
\[ \Delta \text{Deviance} = 0.10131 \]
Hitung Derajat Bebas
Derajat bebas (df) dihitung dengan cara mengurangi derajat bebas model homogenous dengan model saturated:
\[ \text{df} = \text{df}_{\text{model homogenous}} - \text{df}_{\text{model saturated}} = 2 - 0 = 2 \]
Chi-Square Tabel (α = 0.05)
Nilai chi-square untuk derajat bebas 2 pada tingkat signifikansi 5%
dapat dihitung dengan fungsi qchisq:
\[ \chi^2_{0.05, 2} = 5.991465 \]
Keputusan Uji
Keputusan uji dibuat berdasarkan perbandingan antara selisih deviance dan nilai chi-square dari tabel:
\[ \text{Keputusan} = \begin{cases} \text{Terima H0} & \text{jika } \Delta \text{Deviance} \leq \chi^2_{0.05, 2} \\ \text{Tolak H0} & \text{jika } \Delta \text{Deviance} > \chi^2_{0.05, 2} \end{cases} \]
Karena:
\[ \Delta \text{Deviance} = 0.10131 \quad \text{dan} \quad \chi^2_{0.05, 2} = 5.991465 \]
Maka, Terima H0.
Interpretasi Keputusan Uji
Berdasarkan hasil uji, karena \(\Delta \text{Deviance} = 0.10131\) lebih kecil dari \(\chi^2_{0.05, 2} = 5.991465\), kita tidak menolak H0.
Oleh karena itu, kita dapat menyimpulkan bahwa tidak ada bukti yang cukup untuk menunjukkan interaksi tiga arah antara jenis pekerjaan, sikap terhadap pajak tembakau, dan kelompok umur pada tingkat signifikansi 5%.
Model homogenous yang hanya mempertimbangkan interaksi dua arah sudah cukup untuk menjelaskan data, dan tidak perlu menggunakan model saturated dengan interaksi tiga arah.
Uji Model Interaksi Dua Arah (Homogenous vs Conditional on X)
Model Conditional on X
Model log-linear conditional pada X memasukkan efek utama dan interaksi dua arah antara X dengan Y dan X dengan Z, tanpa interaksi antara Y dengan Z maupun interaksi tiga arah.
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} \]
Pengujian Ada Tidaknya Interaksi Antara Y dan Z (Homogenous Model vs Conditional Association on X)
Dalam pengujian ini, kita membandingkan model homogenous dengan model yang mempertimbangkan interaksi antara variabel-variabel dengan mengontrol variabel lain sebagai kondisi (conditional on X).
Hipotesis
Tingkat Signifikansi
Statistik Uji
\[ \Delta \text{Deviance} = \text{Deviance model conditional on X} - \text{Deviance model homogenous} \]
Misalnya, jika hasil deviance untuk kedua model adalah:
\[ \Delta \text{Deviance} = 2.58 - 0.101 = 2.48 \]
\[ db = db_{\text{model conditional on X}} - db_{\text{model homogenous}} = 4 - 2 = 2 \]
Daerah Penolakan
Keputusan
Kesimpulan
Pengujian Selisih Deviance (Conditional on X vs Homogenous)
Deviance of Model
# Model conditional on X
model_conditional_X <- glm(Frequency ~ Job_Type + Tax_Support + Age_Group +
Job_Type*Tax_Support + Job_Type*Age_Group,
family = poisson(link = "log"), data = data)
summary(model_conditional_X)
##
## Call:
## glm(formula = Frequency ~ Job_Type + Tax_Support + Age_Group +
## Job_Type * Tax_Support + Job_Type * Age_Group, family = poisson(link = "log"),
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.6525 0.1260 28.991 < 2e-16 ***
## Job_TypeWorker 0.1795 0.1683 1.067 0.286125
## Tax_SupportSupport 0.7472 0.1280 5.839 5.24e-09 ***
## Age_Group26-40 -0.2877 0.1394 -2.063 0.039105 *
## Age_Group41-60 -0.5390 0.1504 -3.584 0.000339 ***
## Job_TypeWorker:Tax_SupportSupport -0.2772 0.1649 -1.681 0.092840 .
## Job_TypeWorker:Age_Group26-40 0.5108 0.1856 2.752 0.005916 **
## Job_TypeWorker:Age_Group41-60 0.5390 0.1982 2.719 0.006541 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 96.1955 on 11 degrees of freedom
## Residual deviance: 2.5775 on 4 degrees of freedom
## AIC: 87.975
##
## Number of Fisher Scoring iterations: 4
# Selisih deviance antar model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance
Deviance.model
## [1] 2.476187
# Derajat bebas antar model
derajat.bebas <- model_conditional_X$df.residual - model_homogenous$df.residual
derajat.bebas
## [1] 2
# Nilai chi-square tabel dengan alpha = 0.05
chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel
## [1] 5.991465
# Membandingkan selisih deviance dengan nilai chi-square untuk membuat keputusan
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0", "Tolak H0")
Keputusan
## [1] "Terima H0"
Interpretasi:
Karena nilai Deviance.model = 2.48 lebih kecil dari nilai kritis chi-square tabel = 5.99 (dengan df = 2, alpha = 0.05), maka terima H0.
Kesimpulan:
Tidak ada interaksi antara sikap terhadap pajak dan kelompok umur. Dengan kata lain, model yang terbentuk adalah model tanpa parameter interaksi antara Y dan Z.
Uji Model Interaksi Dua Arah (Homogenous vs Conditional on Y)
Model log-linear conditional pada Y memasukkan efek utama dan interaksi dua arah antara X dengan Y dan Y dengan Z, tanpa interaksi antara X dengan Z maupun interaksi tiga arah.
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{jk}^{YZ} \]
Pengujian Ada Tidaknya Interaksi Antara X dan Z (Homogenous Model vs Conditional Association on Y)
Dalam pengujian ini, kita membandingkan model homogenous dengan model yang mempertimbangkan interaksi antara variabel-variabel dengan mengontrol variabel lain sebagai kondisi (conditional on Y).
Hipotesis
H0: Tidak ada interaksi antara Jenis Pekerjaan (Job_Type) dan Kelompok Umur (Age_Group) (model homogenous sudah cukup).
H1: Ada interaksi antara Jenis Pekerjaan (Job_Type) dan Kelompok Umur (Age_Group) (model saturated diperlukan).
Tingkat Signifikansi
\(\alpha = 5\%\)
Statistik Uji
Selisih deviance antara model conditional on Y dan model homogenous:
\[ \Delta \text{Deviance} = \text{Deviance model conditional on Y} - \text{Deviance model homogenous} \]
Misalnya, jika hasil deviance untuk kedua model adalah:
\[ \Delta \text{Deviance} = 10.81 - 0.101 = 10.71 \]
Derajat bebas:
\[ db = db_{\text{model conditional on Y}} - db_{\text{model homogenous}} = 4 - 2 = 2 \]
Daerah Penolakan
Tolak \(H_0\) jika \(\Delta\text{Deviance} > \chi^2_{0.05, 2} = \chi^2_{0.05, 2} = 5.991\)
Keputusan
Karena \(10.71 > 5.991\), maka tolak H0.
Kesimpulan
Dengan taraf nyata 5%, terdapat bukti untuk menolak \(H_0\) atau dapat dikatakan bahwa terdapat interaksi antara jenis pekerjaan dan kelompok umur. Dengan kata lain, model yang terbentuk adalah model dengan parameter interaksi antara Y dan Z.
Pengujian Selisih Deviance (Conditional on Y vs Homogenous)
Deviance of Model
# Model conditional on Y
model_conditional_Y <- glm(Frequency ~ Job_Type + Tax_Support + Age_Group +
Job_Type*Tax_Support + Tax_Support*Age_Group,
family = poisson(link = "log"), data = data)
summary(model_conditional_Y)
##
## Call:
## glm(formula = Frequency ~ Job_Type + Tax_Support + Age_Group +
## Job_Type * Tax_Support + Tax_Support * Age_Group, family = poisson(link = "log"),
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.519e+00 1.344e-01 26.189 < 2e-16 ***
## Job_TypeWorker 5.108e-01 1.333e-01 3.831 0.000128 ***
## Tax_SupportSupport 6.749e-01 1.663e-01 4.058 4.95e-05 ***
## Age_Group26-40 3.222e-16 1.491e-01 0.000 1.000000
## Age_Group41-60 -4.055e-01 1.667e-01 -2.433 0.014983 *
## Job_TypeWorker:Tax_SupportSupport -2.772e-01 1.649e-01 -1.681 0.092840 .
## Tax_SupportSupport:Age_Group26-40 -3.161e-16 1.886e-01 0.000 1.000000
## Tax_SupportSupport:Age_Group41-60 2.624e-01 2.053e-01 1.278 0.201205
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 96.196 on 11 degrees of freedom
## Residual deviance: 10.810 on 4 degrees of freedom
## AIC: 96.207
##
## Number of Fisher Scoring iterations: 4
# Selisih deviance antar model
Deviance.model <- model_conditional_Y$deviance - model_homogenous$deviance
Deviance.model
## [1] 10.70827
# Derajat bebas antar model
derajat.bebas <- model_conditional_Y$df.residual - model_homogenous$df.residual
derajat.bebas
## [1] 2
# Nilai chi-square tabel dengan alpha = 0.05
chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel
## [1] 5.991465
# Membandingkan selisih deviance dengan nilai chi-square untuk membuat keputusan
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0", "Tolak H0")
Keputusan
## [1] "Tolak H0"
Interpretasi:
Karena nilai Deviance.model = 10.71 lebih besar dari nilai kritis chi-square tabel = 5.99 (dengan df = 2, alpha = 0.05), maka tolak H0.
Kesimpulan:
Terdapat interaksi antara jenis pekerjaan dan kelompok umur. Dengan kata lain, model yang terbentuk adalah model dengan parameter interaksi antara X dan Z.
Model Conditional on Z
Model log-linear conditional pada Z memasukkan efek utama dan interaksi dua arah antara X dengan Z dan Y dengan Z, tanpa interaksi antara X dengan Y maupun interaksi tiga arah.
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \]
Pengujian Ada Tidaknya Interaksi Antara X dan Y (Homogenous Model vs Conditional Association on Z)
Dalam pengujian ini, kita membandingkan model homogenous dengan model yang mempertimbangkan interaksi antara variabel-variabel dengan mengontrol variabel lain sebagai kondisi (conditional on Z).
Hipotesis
H0: Tidak ada interaksi antara Jenis Pekerjaan (Job_Type) dan Sikap terhadap pajak (Tax_Support) (model homogenous sudah cukup).
H1: Ada interaksi antara Jenis Pekerjaan (Job_Type) dan Sikap terhadap pajak (Tax_Support) (model saturated diperlukan).
Tingkat Signifikansi
\(\alpha = 5\%\)
Statistik Uji
Selisih deviance antara model conditional on Z dan model homogenous:
\[ \Delta \text{Deviance} = \text{Deviance model conditional on Z} - \text{Deviance model homogenous} \]
Misalnya, jika hasil deviance untuk kedua model adalah:
\[ \Delta \text{Deviance} = 3.32 - 0.101 = 3.22 \]
Derajat bebas:
\[ db = db_{\text{model conditional on Z}} - db_{\text{model homogenous}} = 3 - 2 = 1 \]
Daerah Penolakan
Tolak \(H_0\) jika \(\Delta\text{Deviance} > \chi^2_{0.05, 1} = \chi^2_{0.05, 2} = 3.84\)
Keputusan
Karena \(3.22 < 3.84\), maka terima H0.
Kesimpulan
Dengan taraf nyata 5%, terdapat bukti untuk menerima \(H_0\) atau dapat dikatakan bahwa tidak terdapat interaksi antara jenis pekerjaan dan sikap terhadap pajak. Dengan kata lain, model yang terbentuk adalah model tanpa parameter interaksi antara X dan Y.
Pengujian Selisih Deviance (Conditional on Z vs Homogenous)
Deviance of Model
# Model conditional on Z
model_conditional_Z <- glm(Frequency ~ Job_Type + Tax_Support + Age_Group +
Job_Type*Age_Group + Tax_Support*Age_Group,
family = poisson(link = "log"), data = data)
summary(model_conditional_Z)
##
## Call:
## glm(formula = Frequency ~ Job_Type + Tax_Support + Age_Group +
## Job_Type * Age_Group + Tax_Support * Age_Group, family = poisson(link = "log"),
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.807e+00 1.236e-01 30.797 < 2e-16 ***
## Job_TypeWorker -1.408e-15 1.291e-01 0.000 1.000000
## Tax_SupportSupport 5.108e-01 1.333e-01 3.831 0.000128 ***
## Age_Group26-40 -2.877e-01 1.826e-01 -1.576 0.115095
## Age_Group41-60 -7.108e-01 2.024e-01 -3.512 0.000445 ***
## Job_TypeWorker:Age_Group26-40 5.108e-01 1.856e-01 2.752 0.005916 **
## Job_TypeWorker:Age_Group41-60 5.390e-01 1.982e-01 2.719 0.006541 **
## Tax_SupportSupport:Age_Group26-40 -2.046e-15 1.886e-01 0.000 1.000000
## Tax_SupportSupport:Age_Group41-60 2.624e-01 2.053e-01 1.278 0.201205
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 96.1955 on 11 degrees of freedom
## Residual deviance: 3.3226 on 3 degrees of freedom
## AIC: 90.72
##
## Number of Fisher Scoring iterations: 4
# Selisih deviance antar model
Deviance.model <- model_conditional_Z$deviance - model_homogenous$deviance
Deviance.model
## [1] 3.221263
# Derajat bebas antar model
derajat.bebas <- model_conditional_Z$df.residual - model_homogenous$df.residual
derajat.bebas
## [1] 1
# Nilai chi-square tabel dengan alpha = 0.05
chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel
## [1] 3.841459
# Membandingkan selisih deviance dengan nilai chi-square untuk membuat keputusan
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0", "Tolak H0")
Keputusan
## [1] "Terima H0"
Interpretasi:
Karena nilai Deviance.model = 3.22 lebih kecil dari nilai kritis chi-square tabel = 3.84 (dengan df = 1, alpha = 0.05), maka terima H0.
Kesimpulan:
Tidak terdapat interaksi antara jenis pekerjaan dan sikap terhadap pajak. Dengan kata lain, model yang terbentuk adalah model tanpa parameter interaksi antara X dan Y.
Ringkasan Model Log-Linear
| Model | Parameter | Deviance | Jumlah Parameter | df | AIC |
|---|---|---|---|---|---|
| Saturated | \[ \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} + \lambda_{ijk}^{XYZ} \] | 2.6645e-14 | 12 | 0 | 93.397 |
| Homogenous | \[ \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \] | 0.101 | 10 | 2 | 89.498 |
| Conditional on X | \[ \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} \] | 2.58 | 8 | 4 | 87.975 |
| Conditional on Y | \[ \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{jk}^{YZ} \] | 10.81 | 8 | 4 | 96.207 |
| Conditional on Z | \[ \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \] | 3.32 | 9 | 3 | 90.72 |
Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah
| Interaksi | Pengujian | Δ deviance | Δ df | Chi-square Tabel | Keputusan | Keterangan |
|---|---|---|---|---|---|---|
| XYZ | Saturated vs Homogenous | 0,101 | 2 | 5.991 | Tidak Tolak H0 | tidak ada interaksi |
| YZ | Conditional on X vs Homogenous | 2.48 | 2 | 5.991 | Tidak Tolak H0 | tidak ada interaksi |
| XZ | Conditional on Y vs Homogenous | 10.71 | 2 | 5.991 | Tolak H0 | ada interaksi |
| XY | Conditional on Z vs Homogenous | 3.22 | 1 | 3.841 | Tidak Tolak H0 | tidak ada interaksi |
Kesimpulan Pemilihan Model Terbaik
Dari hasil di atas diketahui bahwa asosiasi yang nyata hanya terdapat antara jenis pekerjaan dan kelompok usia. Sehingga, model terbaik adalah:
\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XZ} \]
Model terbaik adalah model log-linear tanpa interaksi tiga arah dan hanya memuat interaksi dua arah antara jenis pekerjaan dan kelompok usia.
## Model Terbaik
bestmodel <- glm(Frequency ~ Job_Type + Tax_Support + Age_Group + Job_Type*Age_Group,
family = poisson(link = "log"), data = data)
summary(bestmodel)
##
## Call:
## glm(formula = Frequency ~ Job_Type + Tax_Support + Age_Group +
## Job_Type * Age_Group, family = poisson(link = "log"), data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.761e+00 1.049e-01 35.846 < 2e-16 ***
## Job_TypeWorker -1.298e-15 1.291e-01 0.000 1.000000
## Tax_SupportSupport 5.831e-01 8.057e-02 7.237 4.58e-13 ***
## Age_Group26-40 -2.877e-01 1.394e-01 -2.063 0.039105 *
## Age_Group41-60 -5.390e-01 1.504e-01 -3.584 0.000339 ***
## Job_TypeWorker:Age_Group26-40 5.108e-01 1.856e-01 2.752 0.005916 **
## Job_TypeWorker:Age_Group41-60 5.390e-01 1.982e-01 2.719 0.006541 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 96.1955 on 11 degrees of freedom
## Residual deviance: 5.4232 on 5 degrees of freedom
## AIC: 88.82
##
## Number of Fisher Scoring iterations: 4
Berdasarkan perhitungan di atas, didapatkan nilai AIC model terbaik sebesar 88.82. Nilai AIC ini masih lebih besar jika dibandingkan dengan model conditional X. Meskipun begitu, model terbaik yang dipilih tetap model ini karena mempertimbangkan interaksi XZ yang signifikan, sedangkan model conditional X tidak signifikan.
INTERPRETASI KOEFISIEN MODEL TERBAIK
# Interpretasi koefisien model terbaik
data.frame(
koef = bestmodel$coefficients,
exp_koef = exp(bestmodel$coefficients)
)
## koef exp_koef
## (Intercept) 3.760853e+00 42.9850746
## Job_TypeWorker -1.297967e-15 1.0000000
## Tax_SupportSupport 5.831463e-01 1.7916667
## Age_Group26-40 -2.876821e-01 0.7500000
## Age_Group41-60 -5.389965e-01 0.5833333
## Job_TypeWorker:Age_Group26-40 5.108256e-01 1.6666667
## Job_TypeWorker:Age_Group41-60 5.389965e-01 1.7142857
Interpretasi Koefisien Model Terbaik
Berikut adalah interpretasi koefisien berdasarkan output model terbaik yang telah dihitung. Nilai koefisien dan odds ratio memberikan gambaran tentang pengaruh variabel terhadap outcome yang diamati, setelah mengendalikan variabel lain dalam model.
Intercept:
\[ \exp(\lambda_{\text{Intercept}}) =
\exp(3.760853) = 42.99 \]
Tanpa memperhatikan jenis pekerjaan, sikap terhadap pajak, dan kelompok
umur, rata-rata jumlah kasus untuk kategori referensi (misalnya,
Mahasiswa yang tidak mendukung pajak tembakau, dalam kelompok umur 18-25
tahun) adalah 42,99 kali lebih tinggi dibandingkan baseline.
Job_TypeWorker:
\[ \exp(\lambda_{\text{Job\_TypeWorker}}) =
\exp(-1.297967e-15) = 1.00 \]
Tanpa memperhatikan kelompok umur dan sikap terhadap pajak tembakau,
tidak ada perbedaan dalam peluang antara pekerja dan mahasiswa (nilai
odds = 1). Ini menunjukkan bahwa perbedaan jenis pekerjaan (pekerja vs
mahasiswa) tidak memengaruhi hasil yang diamati secara
signifikan.
Tax_SupportSupport:
\[ \exp(\lambda_{\text{Tax\_SupportSupport}})
= \exp(0.5831463) = 1.79 \]
Tanpa memperhatikan jenis pekerjaan dan kelompok umur, peluang seseorang
yang mendukung pajak tembakau adalah 1.79 kali lebih tinggi dibandingkan
yang menolak pajak tembakau.
Age_Group26-40:
\[ \exp(\lambda_{\text{Age\_Group26-40}}) =
\exp(-0.2876821) = 0.75 \]
Tanpa memperhatikan jenis pekerjaan dan sikap terhadap pajak tembakau,
peluang seseorang yang berada di kelompok umur 26-40 tahun adalah 0.75
kali lebih rendah dibandingkan dengan kelompok umur 18-25
tahun.
Age_Group41-60:
\[ \exp(\lambda_{\text{Age\_Group41-60}}) =
\exp(-0.5389965) = 0.58 \]
Tanpa memperhatikan jenis pekerjaan dan sikap terhadap pajak tembakau,
peluang seseorang yang berada di kelompok umur 41-60 tahun adalah 0.58
kali lebih rendah dibandingkan dengan kelompok umur 18-25
tahun.
Job_TypeWorker:Age_Group26-40:
\[
\exp(\lambda_{\text{Job\_TypeWorker:Age\_Group26-40}}) = \exp(0.5108256)
= 1.67 \]
Tanpa memperhatikan sikap terhadap pajak tembakau, odds seseorang yang
bekerja dan berada di kelompok umur 26-40 tahun untuk mendukung pajak
tembakau adalah 1.67 kali lebih tinggi dibandingkan odds yang sama untuk
mahasiswa di kelompok umur yang sama.
Job_TypeWorker:Age_Group41-60:
\[
\exp(\lambda_{\text{Job\_TypeWorker:Age\_Group41-60}}) = \exp(0.5389965)
= 1.71 \]
Tanpa memperhatikan sikap terhadap pajak tembakau, odds seseorang yang
bekerja dan berada di kelompok umur 41-60 tahun untuk mendukung pajak
tembakau adalah 1.71 kali lebih tinggi dibandingkan odds yang sama untuk
mahasiswa di kelompok umur yang sama.
Kesimpulan
Dengan demikian, faktor yang paling berpengaruh adalah sikap terhadap pajak tembakau, diikuti oleh interaksi antara jenis pekerjaan dan kelompok umur, sementara jenis pekerjaan dan kelompok umur tidak berpengaruh signifikan terhadap hasil.
NILAI DUGAAN MODEL TERBAIK
## Fitted Values dari Model Terbaik
data.frame(
Age_Group = z.age,
Job_Type = x.job,
Tax_Support = y.tax,
Frequency = counts,
fitted = bestmodel$fitted.values
)
## Age_Group Job_Type Tax_Support Frequency fitted
## 1 18-25 Student Support 80 77.01493
## 2 18-25 Student Oppose 40 42.98507
## 3 18-25 Worker Support 70 77.01493
## 4 18-25 Worker Oppose 50 42.98507
## 5 26-40 Student Support 60 57.76119
## 6 26-40 Student Oppose 30 32.23881
## 7 26-40 Worker Support 90 96.26866
## 8 26-40 Worker Oppose 60 53.73134
## 9 41-60 Student Support 50 44.92537
## 10 41-60 Student Oppose 20 25.07463
## 11 41-60 Worker Support 80 77.01493
## 12 41-60 Worker Oppose 40 42.98507
Interpretasi:
Kelompok yang Mendukung Pajak Tembakau: Untuk kelompok umur 18-25, baik Student maupun Worker memiliki fitted values yang serupa (77 untuk Support). Ini menunjukkan bahwa mendukung pajak tembakau lebih umum di kalangan kedua kelompok pekerjaan dalam kelompok umur ini.
Kelompok yang Menentang Pajak Tembakau: Student yang menentang pajak tembakau di usia 18-25 memiliki jumlah kejadian yang lebih tinggi dibandingkan dengan Worker yang menentang pajak tembakau di kelompok umur yang sama. Fitted value untuk Student, Oppose di kelompok umur 18-25 adalah 42.99, sementara untuk Worker, Oppose adalah 42.99.
Kelompok Umur yang Lebih Tua: Kelompok umur 41-60 memiliki fitted values yang lebih rendah untuk Student, Oppose (25.07) dibandingkan dengan kelompok umur yang lebih muda (18-25 dan 26-40), yang mengindikasikan bahwa di usia yang lebih tua, jumlah kejadian yang menentang pajak tembakau cenderung lebih rendah.
Efek Pekerjaan: Pekerja (Worker), baik dalam kelompok umur 26-40 dan 41-60, cenderung memiliki fitted values yang lebih tinggi untuk mendukung pajak tembakau dibandingkan dengan Student.
PT Kaca Nusantara adalah perusahaan manufaktur kaca terkemuka di Indonesia yang memproduksi berbagai jenis kaca, termasuk kaca jendela bangunan (float-processed dan non-float), kaca jendela kendaraan, kaca wadah, kaca peralatan makan, dan kaca untuk lampu depan kendaraan. Perusahaan memiliki dataset glass yang berisi 214 sampel kaca dengan komposisi kimia (variabel independen: RI (indeks refraksi), Na (natrium), Mg (magnesium), Al (aluminium), Si (silikon), K (kalium), Ca (kalsium), Ba (barium), Fe (besi)) dan kategori jenis kaca (Type: 1 = jendela bangunan float-processed, 2 = jendela bangunan non-float, 3 = jendela kendaraan float-processed, 5 = wadah, 6 = peralatan makan, 7 = lampu depan kendaraan).
Permasalahan: Perusahaan ingin mengoptimalkan proses pengendalian kualitas dengan mengembangkan sistem klasifikasi otomatis untuk mengidentifikasi jenis kaca berdasarkan komposisi kimianya. Tujuannya adalah memastikan bahwa setiap produk kaca diklasifikasikan dengan benar sebelum dikirim ke pelanggan, menghindari kesalahan produksi, dan meminimalkan biaya akibat pengembalian produk. Selain itu, perusahaan ingin memahami komponen kimia mana yang paling berpengaruh dalam membedakan jenis kaca untuk meningkatkan efisiensi formulasi produksi.
Jawab:
Estimasi Model
library(nnet)
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
##
## Attaching package: 'carData'
## The following object is masked from 'package:vcdExtra':
##
## Burt
##
## Attaching package: 'car'
## The following object is masked from 'package:lawstat':
##
## levene.test
## The following object is masked from 'package:DescTools':
##
## Recode
## The following object is masked from 'package:dplyr':
##
## recode
data <- read.csv(file.choose())
head(data)
## RI Na Mg Al Si K Ca Ba Fe Type
## 1 1.52101 13.64 4.49 1.10 71.78 0.06 8.75 0 0.00 1
## 2 1.51761 13.89 3.60 1.36 72.73 0.48 7.83 0 0.00 1
## 3 1.51618 13.53 3.55 1.54 72.99 0.39 7.78 0 0.00 1
## 4 1.51766 13.21 3.69 1.29 72.61 0.57 8.22 0 0.00 1
## 5 1.51742 13.27 3.62 1.24 73.08 0.55 8.07 0 0.00 1
## 6 1.51596 12.79 3.61 1.62 72.97 0.64 8.07 0 0.26 1
data$Type <- as.factor(data$Type)
model_multinom <- multinom(Type ~ RI + Na + Mg + Al + Si + K + Ca + Ba + Fe, data = data)
## # weights: 66 (50 variable)
## initial value 383.436526
## iter 10 value 257.359885
## iter 20 value 181.634208
## iter 30 value 161.554088
## iter 40 value 157.912577
## iter 50 value 154.889493
## iter 60 value 153.706333
## iter 70 value 153.334999
## iter 80 value 152.219340
## iter 90 value 149.994098
## iter 100 value 149.743843
## final value 149.743843
## stopped after 100 iterations
Interpretasi:
Berdasarkan perhitungan di atas, dapat diketahui bahwa proses 0ptimasi sudah berhasil. Model regresi logistik multinomial berhasil mengoptimalkan parameter, ditunjukkan oleh penurunan nilai log-likelihood dari 383.436526 ke 149.743843 setelah 100 iterasi.
summary(model_multinom)
## Call:
## multinom(formula = Type ~ RI + Na + Mg + Al + Si + K + Ca + Ba +
## Fe, data = data)
##
## Coefficients:
## (Intercept) RI Na Mg Al Si
## 2 114.01139 210.99092 -3.5715880 -6.14888398 -0.0777839 -4.4509190
## 3 46.69565 -61.97027 1.6471464 -0.01788714 2.5121161 0.2207149
## 5 19.54782 14.22700 -0.4893655 -3.69586811 10.1611011 -0.5204113
## 6 -14.59763 -21.52840 10.7663636 -7.48120815 34.9748591 -0.9212133
## 7 -33.83528 22.99089 2.4341715 -5.00880431 6.2849258 -0.1495441
## K Ca Ba Fe
## 2 -3.70543961 -4.6895169 -5.757871 2.2610525
## 3 -0.67459086 0.6082768 -2.208131 1.5301451
## 5 0.62817476 -0.4292740 -3.450644 -0.6424633
## 6 -197.82120395 -4.7069924 -149.906448 -407.9088594
## 7 -0.06454676 -2.2076868 -2.475847 -15.9357312
##
## Std. Errors:
## (Intercept) RI Na Mg Al Si K
## 2 0.12813868 0.32065233 0.5037905 0.7083486 1.1946519 0.1279720 1.7656501
## 3 0.05427163 0.09037508 0.6844281 0.9200780 1.4319207 0.1688232 2.1932557
## 5 0.05937308 0.08790452 0.7494205 1.0489632 2.3181174 0.2053111 2.3195447
## 6 0.07611198 0.12924244 7.2730426 2.0108392 0.7110474 1.8953165 0.1117598
## 7 0.15557299 0.38577177 1.0143998 1.1443317 2.0280315 0.2163546 2.3189653
## Ca Ba Fe
## 2 0.4544060 2.468525 2.050503e+00
## 3 0.5593086 4.500601 3.130438e+00
## 5 0.7324720 2.628137 4.489458e+00
## 6 3.5480224 0.024365 7.124074e-13
## 7 0.9487705 2.840221 1.865328e+01
##
## Residual Deviance: 299.4877
## AIC: 399.4877
Interpretasi Koefisien
Koefisien menunjukkan perubahan log-odds untuk berada di kategori Type tertentu (2, 3, 5, 6, 7) dibandingkan Type 1 per unit peningkatan variabel independen, dengan variabel lain konstan.
(Intercept): 114.01 → Log-odds dasar Type 2 vs. Type 1.
RI: 210.99 → +1 unit RI meningkatkan log-odds Type 2 sebesar 210.99 (efek besar, kemungkinan tidak stabil).
Na: -3.57 → +1 unit Na menurunkan log-odds Type 2 sebesar 3.57.
Mg: -6.15 → +1 unit Mg menurunkan log-odds Type 2 sebesar 6.15.
Al: -0.08 → +1 unit Al menurunkan log-odds Type 2 sebesar 0.08.
Si: -4.45 → +1 unit Si menurunkan log-odds Type 2 sebesar 4.45.
K: -3.71 → +1 unit K menurunkan log-odds Type 2 sebesar 3.71.
Ca: -4.69 → +1 unit Ca menurunkan log-odds Type 2 sebesar 4.69.
Ba: -5.76 → +1 unit Ba menurunkan log-odds Type 2 sebesar 5.76.
Fe: 2.26 → +1 unit Fe meningkatkan log-odds Type 2 sebesar 2.26.
(Intercept): 46.70 → Log-odds dasar Type 3 vs. Type 1.
RI: -61.97 → +1 unit RI menurunkan log-odds Type 3 sebesar 61.97.
Na: 1.65 → +1 unit Na meningkatkan log-odds Type 3 sebesar 1.65.
Mg: -0.02 → +1 unit Mg menurunkan log-odds Type 3 sebesar 0.02.
Al: 2.51 → +1 unit Al meningkatkan log-odds Type 3 sebesar 2.51.
Si: 0.22 → +1 unit Si meningkatkan log-odds Type 3 sebesar 0.22.
K: -0.67 → +1 unit K menurunkan log-odds Type 3 sebesar 0.67.
Ca: 0.61 → +1 unit Ca meningkatkan log-odds Type 3 sebesar 0.61.
Ba: -2.21 → +1 unit Ba menurunkan log-odds Type 3 sebesar 2.21.
Fe: 1.53 → +1 unit Fe meningkatkan log-odds Type 3 sebesar 1.53.
(Intercept): 19.55 → Log-odds dasar Type 5 vs. Type 1.
RI: 14.23 → +1 unit RI meningkatkan log-odds Type 5 sebesar 14.23.
Na: -0.49 → +1 unit Na menurunkan log-odds Type 5 sebesar 0.49.
Mg: -3.70 → +1 unit Mg menurunkan log-odds Type 5 sebesar 3.70.
Al: 10.16 → +1 unit Al meningkatkan log-odds Type 5 sebesar 10.16.
Si: -0.52 → +1 unit Si menurunkan log-odds Type 5 sebesar 0.52.
K: 0.63 → +1 unit K meningkatkan log-odds Type 5 sebesar 0.63.
Ca: -0.43 → +1 unit Ca menurunkan log-odds Type 5 sebesar 0.43.
Ba: -3.45 → +1 unit Ba menurunkan log-odds Type 5 sebesar 3.45.
Fe: -0.64 → +1 unit Fe menurunkan log-odds Type 5 sebesar 0.64.
(Intercept): -14.60 → Log-odds dasar Type 6 vs. Type 1.
RI: -21.53 → +1 unit RI menurunkan log-odds Type 6 sebesar 21.53.
Na: 10.77 → +1 unit Na meningkatkan log-odds Type 6 sebesar 10.77.
Mg: -7.48 → +1 unit Mg menurunkan log-odds Type 6 sebesar 7.48.
Al: 34.97 → +1 unit Al meningkatkan log-odds Type 6 sebesar 34.97 (tidak stabil).
Si: -0.92 → +1 unit Si menurunkan log-odds Type 6 sebesar 0.92.
K: -197.82 → +1 unit K menurunkan log-odds Type 6 sebesar 197.82 (tidak realistis).
Ca: -4.71 → +1 unit Ca menurunkan log-odds Type 6 sebesar 4.71.
Ba: -149.91 → +1 unit Ba menurunkan log-odds Type 6 sebesar 149.91 (tidak realistis).
Fe: -407.91 → +1 unit Fe menurunkan log-odds Type 6 sebesar 407.91 (tidak realistis).
(Intercept): -33.84 → Log-odds dasar Type 7 vs. Type 1.
RI: 22.99 → +1 unit RI meningkatkan log-odds Type 7 sebesar 22.99.
Na: 2.43 → +1 unit Na meningkatkan log-odds Type 7 sebesar 2.43.
Mg: -5.01 → +1 unit Mg menurunkan log-odds Type 7 sebesar 5.01.
Al: 6.28 → +1 unit Al meningkatkan log-odds Type 7 sebesar 6.28.
Si: -0.15 → +1 unit Si menurunkan log-odds Type 7 sebesar 0.15.
K: -0.06 → +1 unit K menurunkan log-odds Type 7 sebesar 0.06.
Ca: -2.21 → +1 unit Ca menurunkan log-odds Type 7 sebesar 2.21.
Ba: -2.48 → +1 unit Ba menurunkan log-odds Type 7 sebesar 2.48.
Fe: -15.94 → +1 unit Fe menurunkan log-odds Type 7 sebesar 15.94.
Residual Deviance: 299.4877 → Mengukur deviasi yang tidak dijelaskan model. Perlu dibandingkan dengan model lain.
AIC: 399.4877 → Menunjukkan kecocokan model dengan penalti kompleksitas. Nilai lebih rendah lebih baik, bandingkan dengan model alternatif.
Nilai P-Value dan Interpretasi
z <- summary(model_multinom)$coefficients / summary(model_multinom)$standard.errors
p_values <- (1 - pnorm(abs(z), 0, 1)) * 2
print(p_values)
## (Intercept) RI Na Mg Al Si K
## 2 0 0 1.346701e-12 0.0000000000 9.480863e-01 0.00000000 0.03584986
## 3 0 0 1.610168e-02 0.9844893895 7.936744e-02 0.19108615 0.75840567
## 5 0 0 5.137616e-01 0.0004261228 1.168725e-05 0.01125292 0.78653091
## 6 0 0 1.387903e-01 0.0001988753 0.000000e+00 0.62693370 0.00000000
## 7 0 0 1.641221e-02 0.0000120293 1.941568e-03 0.48944050 0.97779432
## Ca Ba Fe
## 2 0.0000000 0.01967364 0.2701653
## 3 0.2767932 0.62368802 0.6249862
## 5 0.5578338 0.18919571 0.8862074
## 6 0.1846237 0.00000000 0.0000000
## 7 0.0199710 0.38336700 0.3929317
Interpretasi:
Na, Mg, Al, Si, Ca signifikan untuk beberapa kategori (Type 2, Type 5, Type 7).
K, Ba, Fe signifikan untuk Type 2 dan Type 6, tetapi hasil untuk Type 6 tidak stabil (p-value 0).
RI dan Intercept memiliki p-value 0 di semua kategori, kemungkinan karena ketidakstabilan numerik atau efek kuat.
P-value 0 untuk Al, K, Ba, Fe pada Type 6 menunjukkan masalah estimasi karena jumlah observasi kecil (9 observasi). Kategori ini menyebabkan ketidakstabilan model.
Fe tidak signifikan untuk Type 2, Type 3, Type 5, dan Type 7, menunjukkan kontribusi minimal dalam membedakan kategori ini.
K dan Ba juga sering tidak signifikan, terutama untuk Type 3, Type 5, dan Type 7.
Prediksi dan Validasi
predictions <- predict(model_multinom, data, type = "class")
confusion_matrix <- table(data$Type, predictions)
print(confusion_matrix)
## predictions
## 1 2 3 5 6 7
## 1 52 18 0 0 0 0
## 2 19 54 0 1 0 2
## 3 10 7 0 0 0 0
## 5 0 3 0 9 0 1
## 6 0 0 0 0 9 0
## 7 0 2 0 0 0 27
Interpretasi:
Benar: 52 diprediksi sebagai Type 1.
Salah: 18 diprediksi sebagai Type 2.
Akurasi: 52/70 = 74.29%.
Benar: 54 diprediksi sebagai Type 2.
Salah: 19 sebagai Type 1, 1 sebagai Type 5, 2 sebagai Type 7.
Akurasi: 54/76 = 71.05%.
Benar: 0 diprediksi sebagai Type 3.
Salah: 10 sebagai Type 1, 7 sebagai Type 2.
Akurasi: 0/17 = 0% (model gagal memprediksi Type 3).
Benar: 9 diprediksi sebagai Type 5.
Salah: 3 sebagai Type 2, 1 sebagai Type 7.
Akurasi: 9/13 = 69.23%.
Benar: 9 diprediksi sebagai Type 6.
Salah: 0.
Akurasi: 9/9 = 100% (prediksi sempurna, tetapi jumlah observasi kecil).
Benar: 27 diprediksi sebagai Type 7.
Salah: 2 sebagai Type 2.
Akurasi: 27/29 = 93.10%.
# Hitung akurasi
accuracy <- mean(predictions == data$Type)
cat("Akurasi Model:", accuracy, "\n")
## Akurasi Model: 0.7056075
Interpretasi Akurasi = 0.7056075= 70.56%:
Akurasi dihitung sebagai proporsi prediksi benar (diagonal utama) terhadap total observasi: (52 + 54 + 0 + 9 + 9 + 27) / 214 = 151/214 ≈ 0.7056.
Akurasi 70.56% menunjukkan bahwa model memprediksi dengan benar untuk sekitar 70.56% observasi, tetapi performa ini tergolong sedang dan perlu perbaikan, terutama karena kegagalan memprediksi Type 3.
Kesimpulan
Model tersebut masih tergolong kurang bagus jika akan digunakan di perusahaan besar karena masih banyak kesalahan dalam memprediksi. Selain itu, tingkat keakuratannya juga masih tergolong menengah. Hal tersebut mungkin disebabkan karena kurangnya jumlah observasi pada beberapa tipe kaca, sehingga model tersebut menjadi tidak stabil.
Komponen kimia yang paling berpengaruh dalam membedakan jenis kaca adalah Na, Mg, Al, Si, dan Ca. Komponen-komponen tersebut signifikan untuk beberapa kategori, yaitu tipe 2, tipe 5 dan tipe 7. Beberapa komponen kimia seperti K, Ba, dan Fe juga signifikan untuk tipe 2 dan tipe 6.
Agresti, A. (2010). Analysis of Ordinal Categorical Data, (2nd ed.). Wiley.
Agresti, Alan. (2013). Categorical Data Analysis, 3rd Edition. Wiley.
Agresti, A. (2019). Statistical Methods for the Social Sciences. Pearson.
Friendly, M. (2000). Visualizing Categorical Data. Journal of the American Statistical Association, 95(451), 904-916.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). Chapman and Hall/CRC.
Hosmer, D. W., Lemeshow, S., & Sturdivant, R. X. (2013). Applied Logistic Regression. Wiley.
Long, J. S. (1997). Regression Models for Categorical and Limited Dependent Variables. Sage Publications.
Menard, S. (2002). Applied Logistic Regression Analysis (2nd ed.). Sage Publications.
McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models, 2nd Edition. Chapman and Hall/CRC.
Venables, W. N., & Ripley, B. D. (2002). Modern Applied Statistics with S (4th ed.). New York: Springer.