Segala puji dan syukur penulis panjatkan ke hadirat Tuhan Yang Maha Esa atas berkat, rahmat, dan karunia-Nya sehingga penulisan eBook ini yang berjudul “Analisis Data Kategori” dapat diselesaikan dengan baik.

eBook ini disusun sebagai bagian dari proses pembelajaran dalam mata kuliah Analisis Data Kategori yang diasuh oleh Dr. I Gede Nyoman Mindra Jaya, dosen pada Departemen Statistika Universitas Padjadjaran. Sepanjang perkuliahan, penulis memperoleh banyak pemahaman penting, mulai dari konsep dasar hingga teknik lanjutan dalam analisis data kategori. Materi yang dibahas meliputi distribusi probabilitas diskrit, model-model log-linear, regresi logistik (biner, multinomial, dan ordinal), serta berbagai metode inferensi dalam tabel kontingensi.

Tujuan utama dari penyusunan eBook ini adalah untuk menyajikan rangkuman materi secara lebih terstruktur dan aplikatif agar mudah diikuti oleh mahasiswa maupun peneliti yang tertarik pada analisis data kategori. Setiap bagian dilengkapi dengan penjelasan teori, formulasi matematis, ilustrasi kasus, serta penerapan analisis menggunakan software statistik guna mendukung pemahaman praktis.

Penulis menyadari bahwa buku ini masih jauh dari sempurna, baik dari sisi isi maupun penyajiannya. Oleh karena itu, penulis sangat mengharapkan kritik dan saran yang konstruktif demi penyempurnaan di masa mendatang.

Akhir kata, penulis mengucapkan terima kasih yang sebesar-besarnya kepada Dr. I Gede Nyoman Mindra Jaya atas ilmu, arahan, serta motivasi yang sangat berharga. Semoga eBook ini dapat memberikan manfaat bagi siapa pun yang ingin memperdalam wawasan dalam bidang analisis data kategori, baik di dunia akademik maupun dalam konteks penelitian terapan.

Pendahuluan

Dalam dunia statistika, data yang dikumpulkan tidak selalu berbentuk numerik kontinu. Sering kali, peneliti dihadapkan pada data yang bersifat kategorik, yaitu data yang merepresentasikan kelompok atau klasifikasi, seperti jenis kelamin (laki-laki/perempuan), status merokok (ya/tidak), atau tingkat pendidikan (SD/SMP/SMA). Analisis terhadap data semacam ini dikenal dengan istilah analisis data kategorik.

Berbeda dengan data numerik yang dapat dianalisis menggunakan metode regresi linier, data kategorik memerlukan pendekatan khusus yang dapat menangani struktur data yang terbagi ke dalam kategori. Salah satu teknik yang paling umum digunakan adalah regresi logistik, yang memungkinkan pemodelan hubungan antara satu atau lebih variabel bebas dengan variabel dependen yang bersifat kategorik, khususnya biner (misalnya: sakit/sehat, lulus/tidak lulus).

Tujuan dari analisis data kategorik adalah untuk mengidentifikasi pola, hubungan, dan pengaruh antar variabel, serta untuk membuat prediksi berdasarkan model statistik yang dibangun. Teknik seperti uji chi-square, odds ratio, dan regresi logistik sering digunakan dalam proses ini.

Tujuan Analisis Data Kategorik

Mengidentifikasi Hubungan antara Variabel Kategorik

Analisis data kategorik bertujuan untuk mengetahui apakah terdapat hubungan atau asosiasi yang signifikan antara dua atau lebih variabel kategorik. Contohnya, apakah terdapat hubungan antara jenis kelamin dan preferensi belanja online, atau antara status merokok dan kejadian penyakit jantung. Teknik yang umum digunakan untuk tujuan ini antara lain uji chi-square dan analisis tabulasi silang (cross-tabulation).

Menganalisis Pengaruh Variabel Bebas terhadap Variabel Kategorik

Dalam konteks regresi logistik, analisis data kategorik digunakan untuk mengevaluasi sejauh mana variabel-variabel prediktor (bebas) seperti usia, jenis kelamin, atau tingkat pendidikan memengaruhi probabilitas terjadinya suatu peristiwa biner, seperti “lulus/tidak lulus” atau “terinfeksi/tidak terinfeksi”. Hal ini penting untuk memahami faktor-faktor penyebab dan untuk membangun strategi intervensi atau pencegahan.

Membangun Model Prediktif

Salah satu tujuan penting dari analisis data kategorik adalah membangun model yang dapat digunakan untuk memprediksi peluang atau probabilitas suatu kejadian berdasarkan variabel input. Misalnya, model regresi logistik dapat digunakan untuk memperkirakan kemungkinan seseorang membeli produk berdasarkan usia dan status pekerjaan. Model ini sangat bermanfaat dalam pengambilan keputusan berbasis data (data-driven decision making).

Menyederhanakan Interpretasi Melalui Pengelompokan Kategori

Data kategorik sering kali melibatkan sejumlah besar kategori atau kelas. Analisis ini juga bertujuan untuk menyederhanakan dan mengelompokkan kategori-kategori tersebut agar lebih mudah dianalisis dan diinterpretasikan, tanpa kehilangan makna penting dari data. Ini penting dalam penyusunan laporan dan komunikasi hasil kepada pihak non-statistik.

Definisi dan Ruang Lingkup Analisis Data Kategorik

Analisis data kategorik merupakan cabang dari statistik yang berfokus pada pengolahan dan interpretasi data yang bersifat kategorik, yaitu data yang diklasifikasikan ke dalam kelompok atau kategori tertentu, bukan berupa angka kontinu

Nominal vs Ordinal

Nominal : Kategori yang tidak memiliki urutan atau tingkatan (contoh : jenis kelamin, wilayah, dll)
Ordinal : Kategori yang memiliki urutan atau tingkatan, tetapi selisih antar tingkat tidak diketahui atau tidak sama (contoh : Tingkat Kepuasan, Pendidikan,dll)

Biner vs Multi Kategori

Biner : Hanya memiliki 2 kategori (contoh : Ya/Tidak)
Multikategori : Memiliki lebih dari 2 kategori. (Contoh : Pegawai, Wirausaha, Pelajar, dll)

Perbedaan dengan Data Kuantitatif

Data kategorik mengelompokkan objek ke dalam kategori (tanpa atau dengan urutan), sedangkan data kuantitatif menyatakan besaran numerik yang dapat dihitung atau diukur.

Manfaat Analisis Data Kategori dalam Berbagai Bidang

Kesehatan

Mengkaji hubungan antara penyakit dan faktor risiko (misal: merokok vs tidak, jenis kelamin).

Menganalisis hasil diagnosis (positif/negatif).

Pemasaran dan Bisnis

Mengelompokkan preferensi pelanggan (produk disukai/tidak disukai).

Analisis segmentasi pasar berdasarkan kategori demografis.

Sosial dan Politik

Mengolah data survei opini publik (setuju/tidak setuju).
Menganalisis hasil pemilu berdasarkan kategori wilayah atau kelompok usia.

Pendidikan

Menilai hubungan antara metode pembelajaran dan hasil belajar (lulus/tidak lulus).

Klasifikasi tingkat pendidikan siswa.

Metode

Tabel Kontingensi dan Uji Chi-Square

Tabel kontingensi adalah tabel yang menunjukkan distribusi frekuensi dua (atau lebih) variabel kategorik. Uji Chi-Square (χ²) digunakan untuk mengetahui apakah terdapat hubungan yang signifikan antara dua variabel kategorik. Uji ini membandingkan frekuensi yang diamati dengan frekuensi yang diharapkan jika kedua variabel tidak berhubungan (independen).

Regresi Logistik

Regresi logistik digunakan ketika variabel respons (Y) bersifat kategorik, khususnya biner (misal: ya/tidak, sukses/gagal). Model ini memprediksi peluang (probabilitas) suatu kejadian berdasarkan satu atau lebih variabel prediktor (bisa kategorik atau numerik).

Analisis Correspondence (CA)

Analisis Correspondence (CA) adalah teknik eksploratif grafis untuk menganalisis hubungan antara dua (atau lebih) variabel kategorik dalam bentuk tabel kontingensi.

Decision Tree dan Random Forest

Decision Tree

Metode klasifikasi yang membagi data ke dalam cabang-cabang berdasarkan nilai dari variabel prediktor, untuk memprediksi kategori target. Hasil akhirnya menyerupai pohon dengan simpul keputusan.
Random Forest

Gabungan dari banyak pohon keputusan (decision tree) yang dibangun secara acak. Model ini lebih stabil dan akurat dibanding satu decision tree karena mengurangi overfitting.

Distribusi Probabilitas dalam Data Kategori

Bernoulli

Distribusi Bernoulli menggambarkan kejadian dengan dua kemungkinan hasil: sukses (1) atau gagal (0), dengan probabilitas sukses sebesar \(p\) dan gagal sebesar 1−\(p\). Binomial memiliki fungsi probabilitas sebagai berikut :

\[ P(X=x) =p^x(1-p)^{1-x} \]

dengan,

X = Variabel acak Biner (0 atau 1)

Implementasi dalam R :

set.seed(123)
bernoulli_sample <- rbinom(n = 10, size = 1, prob = 0.5) # 10 percobaan Bernoulli
bernoulli_sample

##  [1] 0 1 0 1 1 0 1 1 1 0

Contoh Soal dalam R :

# Probabilitas sukses 
p <- 0.75

# Menghitung probabilitas untuk sukses (X = 1) dan gagal (X = 0)
prob_success <- p
prob_failure <- 1 - p

# Menampilkan hasil
cat("P(X = 1) =", prob_success, "\n")

## P(X = 1) = 0.75

cat("P(X = 0) =", prob_failure, "\n")

## P(X = 0) = 0.25

Binomial

Distribusi Binomial adalah jumlah dari beberapa percobaan Bernoulli yang independen dan identik, sebanyak \(n\) kali, dengan probabilitas sukses tetap \(p\). Distribusi Binomial memiliki fungsi probabilitas sebagai berikut :

\[ P(X=k) =\binom{n}{k}p^k(1-p)^n-k \]

dengan,

X = Jumlah Keberhasilan

n = Jumlah Percobaan

k = Jumlah Berhasil

Implementasi dalam R :

binomial_sample <- rbinom(n = 10, size = 5, prob = 0.5) 
binomial_sample

##  [1] 4 2 3 3 1 4 2 1 2 4

Contoh Soal dalam R :

# Probabilitas sukses
p <- 0.25

# Jumlah percobaan
n <- 10

# Menghitung probabilitas mendapatkan 3 sukses dalam 10 percobaan
prob_3_success <- dbinom(3, size = n, prob = p)
cat("P(X = 3) =", prob_3_success, "\n")

## P(X = 3) = 0.2502823

Didapatkan bahwa probabilitas mendapatkan 3 sukses dari 10 percobaan dengan p = 0.25 adalah 0.2502823

Multinomial

Distribusi Multinomial adalah perluasan dari distribusi binomial untuk lebih dari dua kategori. Menjelaskan hasil dari \(n\) percobaan, masing-masing dengan probabilitas jatuh ke dalam salah satu dari \(k\) kategori. Multinomial memiliki fungsi probabilitas berikut :

\[ P(X1= x1,...,Xk =xk) =\frac{n!}{x1!x2!...x3!}p1^{x1}p2^{x2}...p3^{x3} \] dengan,

Xi = Frekuensi kemunculan kategori ke-i

n = jumlah total percobaan

xi = Jumlah kejadian kategori ke-i

p1 = probabilitas kategori ke-i

Implementasi dalam R :

set.seed(123)
multinomial_sample <- rmultinom(n = 1, size = 8, prob = c(0.3, 0.5, 0.2))
multinomial_sample

##      [,1]
## [1,]    2
## [2,]    3
## [3,]    3

Contoh Soal dalam R :

# Probabilitas untuk tiap kategori
p <- c(0.4, 0.3, 0.3)

# Jumlah percobaan 
n <- 10

# Menghitung probabilitas distribusi multinomial
multinom_prob <- dmultinom(x = c(2, 3, 5), prob = p)
cat("Probabilitas untuk distribusi multinomial:", multinom_prob, "\n")

## Probabilitas untuk distribusi multinomial: 0.02645395

Probabilitas untuk 2 angka pertama, 3 angka kedua, dan 5 angka pertama adalah 0.026

Poisson

Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam suatu interval waktu atau ruang, ketika kejadian tersebut terjadi secara acak dan independen dengan rata-rata kejadian λ per interval. Poisson memiliki fungsi probabilitas sebagai berikut :

\[ P(X=k)=\frac{e^{-\lambda}\lambda^{k}}{k!} \] dengan,

X : Jumlah kejadian dalam interval tertentu

\(\lambda\) = Rata-rata kejadian dalam interval tertentu

k = Jumlah kejadian yang diamati

Implementasi dalam R :

set.seed(123)
poisson_sample <- rpois(7, lambda = 2) 
poisson_sample

## [1] 1 3 2 4 4 0 2

Contoh Soal dalam R :

# Rata-rata jumlah kejadian
lambda <- 3

# Menghitung probabilitas terjadi 4 kejadian
prob_4_occurrences <- dpois(4, lambda)
cat("P(X = 4) =", prob_4_occurrences, "\n")

## P(X = 4) = 0.1680314

Didapatkan bahwa Probabilitas untuk 4 kejadian dengan lambda = 3 adalah 0.168

Kesimpulan

Distribusi Bernoulli digunakan untuk data kategorikal dengan dua kategori (misalnya sukses dan gagal).
Distribusi Multinomial digunakan untuk data kategorikal dengan lebih dari dua kategori (misalnya hasil lemparan dadu).
Distribusi Binomial digunakan untuk data yang menunjukkan jumlah sukses dalam percobaan berulang.
Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam interval waktu atau ruang tertentu

Desain Sampling

Desain sampling adalah cara merancang pengambilan sampel dalam penelitian, yang krusial untuk validitas hasil. Dalam analisis data kategori, desain sampling dibedakan menjadi dua jenis utama:

Prospective Sampling
Retrospective Sampling

Prospective Sampling

Metode ini melibatkan pengamatan ke depan dalam waktu, artinya subjek diamati seiring waktu setelah mereka dipilih. Umumnya digunakan dalam studi kausal.

Eksperimen

Subjek dibagi secara acak ke dalam kelompok perlakuan dan kontrol. Terdapat beberapa eksperimen yang umum digunakan, yaitu :

Simple Random Sampling (SRS): peluang sama bagi tiap individu.

Stratified Random Sampling: populasi dibagi dalam strata, lalu dipilih acak dari tiap strata.

Cluster Sampling: populasi dibagi menjadi kelompok (cluster), kemudian beberapa cluster dipilih secara acak.

Cohort Study

Kelompok individu dengan karakteristik tertentu diikuti dalam jangka waktu tertentu. Metode yang umum digunakan :

Census Sampling : seluruh populasi diikutsertakan

Systematic Sampling: sampel diambil berdasarkan interval tetap dari daftar.

Matched Sampling : setiap subjek dipasangkan dengan subjek lain berdasarkan variabel tertentu.

Retrospective Sampling

Dalam pendekatan ini, data dikumpulkan dari peristiwa yang telah terjadi. Biasanya digunakan untuk mempelajari hubungan sebab-akibat setelah kejadian berlangsung.

Study Case-Control

Membandingkan individu dengan kondisi tertentu (kasus) dengan yang tidak memiliki kondisi tersebut (kontrol). Metode yang umum digunakan :

Purposive Sampling: berdasarkan karakteristik spesifik.

Snowball Sampling: partisipan awal merekomendasikan partisipan lain.

Incidence Density Sampling: mempertimbangkan waktu munculnya kasus saat memilih kontrol.

Cohort Study Retrospective

Data historis digunakan untuk mengelompokkan subjek berdasarkan paparan dan menganalisis hasil. Metode yang umum digunakan :

Convenience Sampling: berdasarkan data yang tersedia.

Quota Sampling: mencerminkan proporsi populasi tertentu.

Case-Based Sampling : berdasarkan karakteristik kasus yang sudah ada.

Tabel Perbandingan Desain Sampling

Jenis Studi	Pendekatan	Metode Sampling	Keuntungan	Kelemahan
Eksperimen	Prospecitve	SRS, Stratified, Cluster	Kontrol tinggi, analisis kausal	Biaya mahal, isu etika
Cohort Study	Prospecitve	Census, Systematic, Matched	Observasi jangka panjang	Butuh waktu lama, risiko partisipan drop out
Case-Control Study	Retrospecitve	Purposive, Snowball, Incidence Density	Efisien untuk kondisi langka, cepat	Rentan bias recall, kontrol variabel terbatas
Cohort Study Retrospetive	Retrospecitve	Convenience, Quota, Case-Basedundefined	Gunakan data historis, biaya rendah	Ketergantungan pada kualitas data historis

Tabel Kontingensi 2x2

Tabel kontingensi 2×2 adalah tabel yang menyajikan hubungan antara dua variabel kategorik yang masing-masing memiliki dua kategori. Tabel ini sering digunakan dalam penelitian statistik, khususnya epidemiologi, untuk melihat asosiasi antara paparan (exposure) dan kejadian (outcome).

Studi Kasus

Dalam dunia pekerjaan, status pekerjaan dipengaruhi oleh beberapa faktor. Salah satu faktornya ialah Jenis Kelamin. Pada umumnya, laki-laki memiliki peluang kerja yang lebih tinggi dibandingkan perempuan. Oleh karena itu Peneliti ingin melakukan penelitian apakah benar Jenis Kelamin memengaruhi Status Pekerjaan Angkatan Kerja pada Jawa Barat dengan mengambil sampel dari angkatan kerja Kota Bandung?

Status Pekerjaan Angkatan Kerja di Kota Bandung menurut Jenis Kelamin pada tahun 2023
	Bekerja	Tidak Bekerja	Total
Laki-Laki	742536	76725	819261
Perempuan	459949	39705	499654
Total	1202485	116430	1318915

Sumber Data : BPS Kota Bandung. KOTA BANDUNG DALAM ANGKA Bandung Municipality in Figures 2024. vol. 44, Bandung, BPS Kota Bandung, 2024,

# Membuat tabel kontingensi 2 arah
data <- array(
  c(742536, 459949, 76725, 39705),
  dim = c(2, 2),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja")
  )
)


# Menampilkan tabel kontingensi
data

##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki  742536         76725
##     Perempuan  459949         39705

Distribusi Peluang dalam Tabel Kontingensi 2 x 2

Peluang Bersama

Peluang bersama dihitung sebagai: \[ P(Ai, Bj) = \frac{nij}{n} \] Perhitungan menggunakan R :

n = sum(data)
P_Bersama = data / n
P_Bersama

##              Status_Pekerjaan
## Jenis_Kelamin   Bekerja Tidak Bekerja
##     Laki-laki 0.5629900    0.05817282
##     Perempuan 0.3487329    0.03010429

Dari hal ini terlihat bahwa kebanyakan angkatan kerja di Bandung didominasi oleh Laki-Laki yang bekerja yaitu sebesar 56.3%. Sedangkan persentase terkecil didapatkan oleh Perepuan yang tidak bekerja sebesar 3%.

Peluang Marginal

Peluang marginal dihitung sebagai: \[ P(Ai) = \frac{ni.}{n}, \quad P(Bj) = \frac{n.j}{n} \] Perhitungan menggunakan R :

P_Marginal_JenisKelamin = rowSums(data)/n
P_Marginal_StatusPekerjaan = colSums(data)/n
P_Marginal_JenisKelamin; P_Marginal_StatusPekerjaan

## Laki-laki Perempuan 
## 0.6211628 0.3788372

##       Bekerja Tidak Bekerja 
##    0.91172289    0.08827711

Dari hal ini bisa terlihat bahwa Angkatan Kerja didominasi oleh Laki-Laki sebesar 62.1%. Terlihat juga bahwa angkatan kerja didominasi oleh Pekerja sebesar 91.1 %.

Peluang Bersyarat

\[ P(Bj | Ai) = \frac{P(Ai, Bj)}{P(Ai)} = \frac{nij}{ni.} \]

Perhitungan Manual :

Hitung Peluang Bersama

\[ P(Laki-Laki, Bekerja) = \frac{742536}{1318915} \]

\[ P(Laki-Laki,Tidak\ Bekerja) = \frac{76725}{1318915} \]

\[ P(Perempuan, Bekerja) = \frac{459949} {1318915} \]

\[ P(Perempuan, Tidak\ Bekerja) = \frac{39705}{1318915} \]
Hitung Peluang Marginal

\[ P(Laki-Laki)=\frac{819261}{1318915} \]

\[ P(Perempuan)=\frac{499654}{1318915} \]

\[ P(Bekerja)=\frac{1202485}{1318915} \]

\[ P(Tidak\ Bekerja) = \frac{116430}{1318915} \]
Hitung Peluang Bersyarat

\[ P(Bekerja|Laki-Laki) = \frac{742536}{819261} \]

\[ P(Kanker|Tidak\ Perokok) = \frac{459949}{499654} \]

Perhitungan menggunakan R :

P_Conditional = data / rowSums(data) 
P_Conditional

##              Status_Pekerjaan
## Jenis_Kelamin   Bekerja Tidak Bekerja
##     Laki-laki 0.9063485    0.09365147
##     Perempuan 0.9205350    0.07946499

Dari Hal ini terlihat bahwa Peluang Perempuan yang bekerja (92%) lebih tinggi dibandingkan dengan Laki-Laki yang Bekerja (90.6%).

Ukuran Asosiasi

Ukuran asosiasi digunakan untuk mengukur seberapa kuat hubungan antara dua variabel kategori, terutama paparan dan kejadian. Tiga ukuran yang umum digunakan:

Perbedaan Peluang (Risk Difference, RD)

RD adalah selisih antara risiko kejadian pada kelompok terpapar dan tidak terpapar.

\[ RD = {P(E|X=1)} - {P(E|X=0)} \]

Interpretasi:

RD > 0: risiko lebih besar pada kelompok 1.
RD < 0: risiko lebih kecil pada kelompok 1.
RD = 0: tidak ada perbedaan risiko antara kedua kelompok.

Perhitungan menggunakan R :

RD = function(n11, n12, n21, n22) {
(n11 / (n11 + n12)) - (n21 / (n21 + n22))
}
RD(742536, 76725, 459949, 39705)

## [1] -0.01418648

Dari hal diatas dapat dikatakan bahwa perempuan memiliki risiko 1.4% lebih tinggi dalam bekerja dibandingkan laki-laki.

Risiko Relatif (Relative Risk, RR)

RR membandingkan proporsi kejadian antara dua kelompok.

\[ RR = \frac{P(E|X=1)}{P(E|X=0)} \]

Interpretasi:

RR > 1: risiko lebih tinggi pada kelompok 1.
RR < 1: risiko lebih rendah pada kelompok 1.
RR = 1: tidak ada perbedaan risiko antar kelompok.

Perhitungan Menggunakan R :

RR = function(n11, n12, n21, n22) {
(n11 / (n11 + n12)) / (n21 / (n21 + n22))
}
RR(742536, 76725, 459949, 39705)

## [1] 0.9845889

Dari hal diatas dapat dikatakan bahwa laki-laki memliki risiko 0.98 lebih kecil mendapatkan pekerjaan dibandingkan perempuan.

Odds Ratio (OR)

OR membandingkan peluang terjadinya suatu kejadian antara dua kelompok. Digunakan terutama dalam studi kasus-kontrol.

\[ OR = \frac{a \times d}{b \times c} \]

Interpretasi:

OR > 1: peluang kejadian lebih besar pada kelompok 1.
OR < 1: peluang lebih kecil pada kelompok 1.
OR = 1: tidak ada perbedaan peluang.

Perhitungan menggunakan R :

OR = function(n11, n12, n21, n22) {
(n11*n22)/(n12*n21)
  }
OR(742536, 76725, 459949, 39705)

## [1] 0.8354417

Dari hal diatas dapat disimpulkan bahwa Laki-Laki memiliki 0.83 lebih kecil odds mendapatkan pekerjaan dibandingkan Perempuan.

Ukuran	Fokus	Studi yg cocok	Interpretasi Utama
RD	Selisih Risiko	Cohort, Eksperimen	Tambahan risiko atau perlindungan akibat paparan
RR	Proporsi Risiko	Cohort, Eksperimen	Seberapa kali risiko lebih besar pada kelompok terpapar
OR	Perbandingan Odds	Case-Control	Estimasi risiko relatif dalam studi non-prospektif

Inferensi Tabel Kontingensi 2x2

Inferensi statistik dilakukan untuk mengambil kesimpulan terhadap populasi berdasarkan data sampel. Dalam konteks tabel kontingensi 2×2, inferensi terbagi menjadi estimasi, pengujian hipotesis, dan analisis residual.

Estimasi

Estimasi merupakan proses penarikan kesimpulan tentang parameter populasi berdasarkan data sampel.

Estimasi Titik

Estimasi titik memberikan nilai tunggal sebagai taksiran parameter populasi. Misalnya:

\[ p = \frac{x}{n} \]

dengan, p = estimasi titik proporsi, x = jumlah individu dalam kategori tertenu, dan n = total jumlah individu dalam sampel.

Estimasi Interval

Estimasi interval memberikan rentang nilai untuk parameter dengan tingkat kepercayaan tertentu.

\[ p\pm Z\alpha/2 \sqrt{\frac{p(1-p)}{n}} \]

dengan,

\(Z\alpha/2\)= nilai dari distribusi normal standar untuk tingkat kepercayaan tertentu p = estimasi titik proporsi

n = ukuran sampel

Uji Hipotesis

Uji Proporsi

Uji proporsi dua kelompok bertujuan untuk mengetahui apakah terdapat perbedaan proporsi kejadian antara dua kelompok. Uji ini dilakukan dengan menggunakan uji statistik z dua proporsi.

Estimasi proporsi kelompok sebagai berikut :

\[ p1 = \frac{n11}{n1.} \]

\[ p2 = \frac{n21}{n2.} \]

Estimasi proporsi gabungan :

\[ p = \frac{n11+n21}{n1.+n2.} \]

Statistik Uji :

\[ Z = \frac{p1-p2}{\sqrt(p(1-p)(\frac{1}{n1.}+\frac{1}{n2.})} \]

Hipotesis

H₀ : Tidak ada perbedaan proporsi antara Laki-Laki dan Perempuan (p₁ = p₂)

H₁ : Tidak ada perbedaan proporsi antara Laki-Laki dan Perempuan (p₁ \(\ne\) p₂)
Taraf Signifikansi

\(\alpha\) = 5 %
Statistik Uji

# Membuat tabel kontingensi 2 arah
data <- array(
  c(742536, 459949, 76725, 39705),
  dim = c(2, 2),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja")
  )
)

prop_test <- prop.test(x = c(data[1,1], data[2,1]),
n = c(sum(data[1,]), sum(data[2,])))
print(prop_test)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(data[1, 1], data[2, 1]) out of c(sum(data[1, ]), sum(data[2, ]))
## X-squared = 775.92, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.01516809 -0.01320487
## sample estimates:
##    prop 1    prop 2 
## 0.9063485 0.9205350

Kriteria Uji

Tolak H₀ jika p-value < \(\alpha\)
Keputusan

Tolak H₀ karena p-value < \(\alpha\) (2.2e-16 < 0.05)
Kesimpulan

Dapat disimpulkan bahwa Tidak ada perbedaan proporsi antara Laki-Laki dan Perempuan (p₁ \(\ne\) p₂) setelah dilakukan Uji Proporsi dengan \(\alpha\) = 5%.

Uji Asosiasi

Mengukur kekuatan hubungan antara dua variabel kategori.

RD

Standar Error untuk RD diperoleh sebagai berikut :

\[ SE(RD) = \sqrt{(\frac{p1(1-p1)}{n1.})+(\frac{p2(1-p2)}{n2.})} \]

Statistik Uji :

\[ Z = \frac{RD}{SE(RD)} \]

Perhitungan Menggunakan R :

n11 = 742536
n12 = 76725
n21 = 459949
n22 = 39705
n1. = n11 + n12
n2. = n21 + n22

# Risk Difference
p1<-(n11/n1.)
p2<-(n21/n2.)
rd <- p1 - p2
se_rd <- sqrt((p1 * (1 - p1) / n1.) + p2*((1 - p2) / n2.))
z_rd <- rd / se_rd
list(RD =rd ,SE_RD = se_rd, Z_RD = z_rd)

## $RD
## [1] -0.01418648
## 
## $SE_RD
## [1] 0.0005000086
## 
## $Z_RD
## [1] -28.37248

RR

Standar Error untuk RR diperoleh sebagai berikut :

\[ SE(ln RR) = \sqrt{(\frac{1}{n11})-(\frac{1}{n1.})+(\frac{1}{n21})-(\frac{1}{n2.})} \]

Statistik Uji :

\[ Z = \frac{ln RR}{SE(ln RR)} \]

Perhitungan Menggunakan R :

# Relative Risk
rr <- (n11/n1.) / (n21/n2.)
se_ln_rr <- sqrt((1/n11) - (1/n1.) + (1/n21) - (1/n2.))
z_rr <- log(rr) / se_ln_rr
list(RR =rr ,SE_RR = se_ln_rr, Z_RR = z_rr)

## $RR
## [1] 0.9845889
## 
## $SE_RR
## [1] 0.000546711
## 
## $Z_RR
## [1] -28.40827

OR

Standar Error untuk OR diperoleh sebagai berikut :

\[ SE(ln OR) = \sqrt{(\frac{1}{n11})+(\frac{1}{n12})+(\frac{1}{n21})+(\frac{1}{n22})} \]

Statistik Uji :

\[ Z = \frac{ln OR}{SE(ln OR)} \]

Perhitungan menggunakan R :

# Odds Ratio
or <- (n11 * n22) / (n12 * n21)
se_ln_or <- sqrt((1/n11) + (1/n12) + (1/n21) + (1/n22))
z_or <- log(or) / se_ln_or
list(OR = or, SE_Ln_OR = se_ln_or, Z_OR = z_or)

## $OR
## [1] 0.8354417
## 
## $SE_Ln_OR
## [1] 0.006460665
## 
## $Z_OR
## [1] -27.82914

Uji Independensi

Uji independensi digunakan untuk menentukan apakah ada hubungan statistik antara dua variabel kategorikal.

Chi-Square

Uji Chi-Square digunakan untuk menguji apakah ada hubungan antara dua variabel kategorikal

Statistik Uji :

\[ \chi^2 = \sum\frac{(O-E)^2}{E} \]

dengan,

O = Nilai Observasi dalam tabel kontingensi

E = Nilai yang diharapkan dengan

\[Eij = \frac{(Ri \times Cj)}{N}\]

dengan,

Ri = Total Baris ke i

Cj = Total Kolom ke j

N = Total Data

Langkah Manual :

Hitung Frekuensi Harapan untuk masing-masing sel

\[ E_{11}=\frac{819261 \times 1202485}{1318915} \]

\[ E_{21}=\frac{499654\times 1202485}{1318915} \]

\[ E_{12}=\frac{819261 \times 116430}{1318915} \]

\[ E_{22}=\frac{499654\times 116430}{1318915} \]
Hitung Nilai Chi-Square

Perhitungan Menggunakan R :

# Membuat tabel kontingensi 2 arah
data <- array(
  c(742536, 459949, 76725, 39705),
  dim = c(2, 2),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja")
  )
)

# Uji Chi-Square
chisq_test <- chisq.test(data)
print(chisq_test)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data
## X-squared = 775.92, df = 1, p-value < 2.2e-16

Karena p-value < \(\alpha\) (2.2e-16 < 0.05) maka didaptkan bahwa terdapat hubungan Jenis Kelamin dan Status Pekerjaan.

Partisi Chi-Square

Membagi nilai chi-square menjadi bagian-bagian untuk menganalisis efek interaksi atau kontribusi tiap kategori.

Pada kolom Tidak bekerja dapat dipecah menjadi 2 kolom yaitu “Pernah Bekerja” dan “Tidak Pernah Bekerja” seperti berikut :

Status Pekerjaan Angkatan Kerja di Kabupaten Bandung Barat menurut Jenis Kelamin pada tahun 2023
	Bekerja	Pernah Bekerja	Tidak Pernah Bekerja	Total
Laki-Laki	575757	28300	26948	631005
Perempuan	294512	8146	13461	316119
Total	870269	36446	40409	947124

Sumber Data : https://bandungbaratkab.bps.go.id/id/statistics-table/1/Mjc0IzE=/perkembangan-penduduk-kabupaten-bandung-barat-berumur-15-tahun-ke-atas-dan-jenis-kegiatan-selama-seminggu-yang-lalu-menurut-jenis-kelamin--agustus-2023.html

Perhitungan Menggunakan R :

# Membuat tabel kontingensi 2 arah
data1 <- array(
  c(575757
, 294512
, 28300
, 8146
, 26948
, 13461
),
  dim = c(2, 3),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Pernah Bekerja","Tidak Pernah Bekerja")
  )
)


# Menampilkan tabel kontingensi
data1

##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Pernah Bekerja Tidak Pernah Bekerja
##     Laki-laki  575757          28300                26948
##     Perempuan  294512           8146                13461

# Uji Chi-Square
chisq_test1 <- chisq.test(data1)
print(chisq_test1)

## 
##  Pearson's Chi-squared test
## 
## data:  data1
## X-squared = 2077.1, df = 2, p-value < 2.2e-16

Karena p-value <\[\alpha\] (2.2e-16 < 0.05) maka didaptkan bahwa terdapat hubungan Jenis Kelamin dan Status Pekerjaan.

Likelihood Ratio

Uji Likelihood Ratio (G²) adalah alternatif dari uji chi-square yang digunakan untuk menguji hipotesis independensi dalam tabel kontingensi I×J.

Statistik Uji :

\[ G^2 = 2\sum\sum nij \times ln(\frac{nij}{\mu ij}) \]

dengan,

nij = frekuensi observasi

\[ \mu ij = n \times pi. \times p.j = frekuensi ekspektasi \]

Perhitungan menggunakan R :

# Hitung Frekuensi Ekspektasi
data_expected <- chisq.test(data)$expected

# Hitung Statistik G²
G2 <- 2 * sum(data * log(data / data_expected))

# Nilai kritis 
critical_value <- qchisq(0.95, df = 1)

list(G2 = G2, Titik_Kritis = critical_value)

## $G2
## [1] 786.1287
## 
## $Titik_Kritis
## [1] 3.841459

Karena G^2 > \[\chi^2\] (786.128 > 3.841) maka didapatkan bahwa terdapat hubungan antara Jenis Kelamin dan Status Pekerjaan

Exact Fisher

Uji Fisher’s Exact digunakan untuk menguji hubungan antara dua variabel kategorikal dalam tabel kontingensi kecil, dimana asumsi Chi-square tidak berlaku karena ukuran sampel yang kecil.

\[ P = \frac{(a+b)!\times(c+d)!\times(a+c)!\times (b+d)!}{a!\times b!\times c! \times d!} \]

dengan :

a,b,c,d = frrekuensi dalam sel

n = total observasi

Perhitungan Menggunakan R :

# Membuat tabel kontingensi 2 arah
data <- array(
  c(742536, 459949, 76725, 39705),
  dim = c(2, 2),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja")
  )
)

fisher.test(data)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  data
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.8249179 0.8460769
## sample estimates:
## odds ratio 
##  0.8354544

Karena p-value < \[\alpha\] (2.2e-16 < 0.05) maka didaptkan bahwa terdapat hubungan Jenis Kelamin dan Status Pekerjaan.

Analisis Residual

Residual digunakan untuk menilai kontribusi setiap sel dalam tabel kontingensi terhadap hubungan antara variabel. Residual mengukur selisih antara frekuensi observasi (O) dan frekuensi ekspektasi (E) berdasarkan asumsi independensi.

Residual besar (positif/negatif): menunjukkan sel menyumbang signifikan terhadap asosiasi.
Residual mendekati 0: data sesuai dengan asumsi independensi; tidak ada asosiasi yang kuat.

Jenis Residual

Pearson Residual

\(Rij =\frac{Oij-Eij}{\sqrt{Eij}}\)

Dengan,

Oij = nilai observasi pada sel ke i ,j

Eij = nilai ekspektas pada sel ke i,j
Standardized Residual

\(rij =\frac{Oij-Eij}{\sqrt{Eij(1-pi+)(1-p+j)}}\)

Perhitungan Residual Manual

Contoh Soal :

	Kategori A	Kategori B	Total
Grup 1	20	10	30
Grup 2	30	20	50
Total	50	30	80

Hitung Frekuensi yang diharapkan

\[ Eij=\frac{(Total\ Baris) \times(Toal Kolom)}{Total\ Keseluruhan} \] dengan,
- \[ E11 = \frac{30 \times50}{80} \]
  
  \[ E12 = \frac{30 \times30}{80} \]
  
  \[ E21 = \frac{50 \times50}{80} \]
  
  \[ E22 = \frac{50 \times30}{80} \]
Hitung Pearson Residual

\[ e11 = \frac{20-18.75}{\sqrt{18.75}} = 0.29 \]

\[ e12 = \frac{10-11.25}{\sqrt{11.25}} = -0.37 \]

\[ e21 = \frac{30-31.25}{\sqrt{31.25}} = -0.22 \]

\[ e22 = \frac{20-18.75}{\sqrt{18.75}} = 0.29 \]

Perhitungan Menggunakan R :

expected <- chisq.test(data)$expected
pearson_residual = (data - expected) / sqrt(expected)

row_sum <- rowSums(data)
col_sum <- colSums(data)
total_sum <- sum(data)

standardized_residual <- (data - expected) / 
  sqrt(expected * (1 - row_sum / total_sum) %*% t(1 - col_sum / total_sum))


list(Pearson_Residual = pearson_residual, Standarized_Residual = standardized_residual)

## $Pearson_Residual
##              Status_Pekerjaan
## Jenis_Kelamin   Bekerja Tidak Bekerja
##     Laki-laki -5.094564      16.37248
##     Perempuan  6.523543     -20.96481
## 
## $Standarized_Residual
##              Status_Pekerjaan
## Jenis_Kelamin   Bekerja Tidak Bekerja
##     Laki-laki -27.85844      27.85844
##     Perempuan  27.85844     -27.85844

Deteksi Outlier menggunakan Residual

Apa Itu Outlier dalam Tabel Kontingensi?

Outlier adalah sel dengan residual besar, menunjukkan frekuensi yang menyimpang jauh dari ekspektasi (baik jauh lebih tinggi maupun lebih rendah)

Menggunakan Residual untuk Mendeteksi Outlier

Pearson Residual → Jika |e| > 2 → kemungkinan ada outlier.

Standardized Residual → Jika |r| > 3 → indikasi kuat adanya outlier signifikan.

Perhitungan menggunakan R

# Menampilkan hasil
list(
Observed = data,
Expected = expected,
Pearson_Residual = pearson_residual,
Standardized_Residual = standardized_residual
)

## $Observed
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki  742536         76725
##     Perempuan  459949         39705
## 
## $Expected
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki  746939      72321.99
##     Perempuan  455546      44108.01
## 
## $Pearson_Residual
##              Status_Pekerjaan
## Jenis_Kelamin   Bekerja Tidak Bekerja
##     Laki-laki -5.094564      16.37248
##     Perempuan  6.523543     -20.96481
## 
## $Standardized_Residual
##              Status_Pekerjaan
## Jenis_Kelamin   Bekerja Tidak Bekerja
##     Laki-laki -27.85844      27.85844
##     Perempuan  27.85844     -27.85844

Karena Residual dari Pearson > 3 didaparkan bahwa sel-sel tersebut terdapat outlier yang signifikan. Selain itu residual besar menunjukkan hubungan signifikan antar variabel

Tabel Kontingensi 3 Arah

Tabel kontingensi tiga arah digunakan untuk menganalisis tiga variabel kategorik secara simultan, misalnya untuk menilai hubungan antara (X), (Y), dan (Z). Variabel ketiga (Z) disebut sebagai variabel kontrol karena digunakan untuk mengendalikan efek terhadap hubungan antara X dan Y

Studi Kasus

Dalam dunia pekerjaan, status pekerjaan dipengaruhi oleh beberapa faktor. Salah satu faktornya ialah Jenis Kelamin. Pada umumnya, laki-laki memiliki peluang kerja yang lebih tinggi dibandingkan perempuan. Namun pada beberapa kasus, hal ini bisa berbanding terbalik dimana jumlah pekerja Perempuan lebih banyak daripada Laki-Laki. Banyak juga faktor yang memengaruhi hal tersebut, salah satunya wilayah.

Dengan kata lain, wilayah sesorang juga dapat berpengaruh terhadap status pekerjaan. Oleh karena itu Peneliti ingin melihat apakah benar bahwa Jenis Kelamin dan Wilayah memengaruhi Status Pekerjaan seseorang ? Peneliti akan melakukan penelitian di 4 wilayah Bandung Raya, yaitu Kota Bandung, Kabupaten Bandung, Kota Cimahi, dan Kabupaten Bandung Raya.

Status Pekerjaan Angkatan Kerja di Wilayah Bandung Raya menurut Jenis Kelamin pada tahun 2023.
Z (Wilayah)	X (Jenis Kelamin)	Y (Bekerja)	Y (Tidak Bekerja)	Total
Kota Bandung	Laki-Laki	742536	76725	819261
	Perempuan	459949	39705	499654
Kabupaten Bandung	Laki-Laki	1116497	84723	1201220
	Perempuan	639831	37706	677537
Kota Cimahi	Laki-Laki	166488	27298	193786
	Perempuan	115910	5894	121804
Kabupaten Bandung Barat	Laki-Laki	575757	55248	631005
	Perempuan	294512	21607	316119
Total		4111480	348906	4460386

Sumber Data :

BPS Kabupaten Bandung. KABUPATEN BANDUNG DALAM ANGKA Bandung Regency in Figures 2024. vol. 41, Kabupaten Bandung, BPS Kabupaten Bandung, 2024,
BPS Kabupaten Bandung Barat. KABUPATEN BANDUNG BARAT DALAM ANGKA Bandung Barat Regency in Figures 2024. vol. 15, Kabupaten Bandung Barat, BPS Kabupaten Bandung Barat, 2024,
BPS Kota Bandung. KOTA BANDUNG DALAM ANGKA Bandung Municipality in Figures 2024. vol. 44, Bandung, BPS Kota Bandung, 2024,
BPS Kota Cimahi. KOTA CIMAHI DALAM ANGKA Cimahi Municipality in Figures 2024. vol. 2022, Kota Cimahi, BPS Kota Cimahi, 2024,

Implementasi dalam R

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)


# Menampilkan tabel kontingens
data2

## , , Wilayah = Kota Bandung
## 
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki  742536         76725
##     Perempuan  459949         39705
## 
## , , Wilayah = Kabupaten Bandung
## 
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki 1116497         84723
##     Perempuan  639831         37706
## 
## , , Wilayah = Kota Cimahi
## 
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki  166488         27298
##     Perempuan  115910          5894
## 
## , , Wilayah = Kabupaten Bandung Barat
## 
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki  575757         55248
##     Perempuan  294512         21607

Tabel Parsial dan Marginal

Tabel Frekuensi Parsial

Tabel ini menyajikan hubungan antara dua variabel (X dan Y) dalam setiap wilayah Z.

Implementasi dalam R :

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

# Ekstrak tabel parsial berdasarkan wilayah
freq_parsial_kotabandung <- data2[, , "Kota Bandung"]
freq_parsial_kabbandung <- data2[, , "Kabupaten Bandung"]
freq_parsial_kotacimahi <- data2[, , "Kota Cimahi"]
freq_parsial_kabbandungbrt <- data2[, , "Kabupaten Bandung Barat"]
# Tampilkan hasil
list(Kota_Bandung = freq_parsial_kotabandung,Kabupaten_Bandung = freq_parsial_kabbandung, Kota_Cimahi = freq_parsial_kotacimahi,Kabupaten_Bandung_Barat = freq_parsial_kabbandungbrt)

## $Kota_Bandung
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki  742536         76725
##     Perempuan  459949         39705
## 
## $Kabupaten_Bandung
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki 1116497         84723
##     Perempuan  639831         37706
## 
## $Kota_Cimahi
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki  166488         27298
##     Perempuan  115910          5894
## 
## $Kabupaten_Bandung_Barat
##              Status_Pekerjaan
## Jenis_Kelamin Bekerja Tidak Bekerja
##     Laki-laki  575757         55248
##     Perempuan  294512         21607

Tabel Frekuensi Marginal

Frekuensi marginal dihitung dengan menjumlahkan data dari semua kategori variabel gender (X) dan wilayah (Z).

Implementasi dalam R :

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

# Hitung frekuensi marginal
freq_marginal_X <- apply(data2, 1, sum)
freq_marginal_Z <- apply(data2, 3, sum)
# Tampilkan hasil
list(Peluang_Marginal_Gender = freq_marginal_X, Peluang_Marginal_Wilayah = freq_marginal_Z)

## $Peluang_Marginal_Gender
## Laki-laki Perempuan 
##   2845272   1615114 
## 
## $Peluang_Marginal_Wilayah
##            Kota Bandung       Kabupaten Bandung             Kota Cimahi 
##                 1318915                 1878757                  315590 
## Kabupaten Bandung Barat 
##                  947124

Distribusi Peluang

Peluang Bersama

\[ P(Z,X,Y) =\frac{f(Z,X,Y)}{N} \]

Perhitungan menggunakan R :

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

P_ZXY = data2 / sum(data2)
P_ZXY

## , , Wilayah = Kota Bandung
## 
##              Status_Pekerjaan
## Jenis_Kelamin   Bekerja Tidak Bekerja
##     Laki-laki 0.1664735   0.017201426
##     Perempuan 0.1031187   0.008901696
## 
## , , Wilayah = Kabupaten Bandung
## 
##              Status_Pekerjaan
## Jenis_Kelamin   Bekerja Tidak Bekerja
##     Laki-laki 0.2503140   0.018994544
##     Perempuan 0.1434475   0.008453528
## 
## , , Wilayah = Kota Cimahi
## 
##              Status_Pekerjaan
## Jenis_Kelamin    Bekerja Tidak Bekerja
##     Laki-laki 0.03732592   0.006120098
##     Perempuan 0.02598654   0.001321410
## 
## , , Wilayah = Kabupaten Bandung Barat
## 
##              Status_Pekerjaan
## Jenis_Kelamin    Bekerja Tidak Bekerja
##     Laki-laki 0.12908233    0.01238637
##     Perempuan 0.06602837    0.00484420

Peluang Marginal

Diperoleh dengan menjumlahkan semua nilai peluang bersama untuk variabel lain.

Misalnya Peluang Bekerja :

\[ P(Y=Bekerja) = \sum P(Z,X,Y=Bekerja) \]

Perhitungan Menggunakan R :

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

# Hitung probabilitas marginal
marginal_X <- apply(P_ZXY, 1, sum)
marginal_Z <- apply(P_ZXY, 3, sum)
# Tampilkan hasil
list(Peluang_Marginal_Gender = marginal_X, Peluang_Marginal_Wilayah =
marginal_Z)

## $Peluang_Marginal_Gender
## Laki-laki Perempuan 
## 0.6378982 0.3621018 
## 
## $Peluang_Marginal_Wilayah
##            Kota Bandung       Kabupaten Bandung             Kota Cimahi 
##              0.29569526              0.42120951              0.07075397 
## Kabupaten Bandung Barat 
##              0.21234126

Peluang Bersyarat

\[ P(Z|X,Y)=\frac{Z,X,Y}{X,Y} \]

Sebagai contoh, peluang seseorang perempuan diberikan bahwa dia bekerja dan tinggal di Cimahi :

\[ P(Cimahi|X=Perempuan, Y=Bekerja) = \frac{P(Cimahi,Perempuan, Bekerja)}{P(Perempuan, Bekerja)} \]

Perhitungan menggunakan R :

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

# Ambil data untuk Kota Cimahi
data_cimahi <- data2[ , , "Kota Cimahi"]

# Jumlah orang yang BEKERJA di Cimahi (laki-laki & perempuan)
total_bekerja_cimahi <- sum(data_cimahi[ , "Bekerja"])

# Jumlah PEREMPUAN yang BEKERJA di Cimahi
perempuan_bekerja_cimahi <- data_cimahi["Perempuan", "Bekerja"]

# Hitung peluang bersyarat
P_perempuan_given_bekerja_cimahi <- perempuan_bekerja_cimahi / total_bekerja_cimahi

# Tampilkan hasil
P_perempuan_given_bekerja_cimahi

## [1] 0.4104491

Tabel Peluang Bersyarat

Peluang bersyarat dihitung berdasarkan jumlah total dalam setiap kategori.

Perhitungan menggunakan R :

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

ftable(data2)

##                                Wilayah Kota Bandung Kabupaten Bandung Kota Cimahi Kabupaten Bandung Barat
## Jenis_Kelamin Status_Pekerjaan                                                                           
## Laki-laki     Bekerja                        742536           1116497      166488                  575757
##               Tidak Bekerja                   76725             84723       27298                   55248
## Perempuan     Bekerja                        459949            639831      115910                  294512
##               Tidak Bekerja                   39705             37706        5894                   21607

Ukuran Asosiasi

Perbedaan Peluang (Risk Difference, RD)

\[ RD = P(Y|X1,Z) - P(Y|X2,Z) \]

Perhitungan menggunakan R :

#Z = Kota Bandung
p1 <- freq_parsial_kotabandung[1, 1] / sum(freq_parsial_kotabandung[1, ])
p2 <- freq_parsial_kotabandung[2, 1] / sum(freq_parsial_kotabandung[2, ])
RDKotaBandung <- p1 - p2

#Z = Kabupaten Bandung
p1 <- freq_parsial_kabbandung[1, 1] / sum(freq_parsial_kabbandung[1, ])
p2 <- freq_parsial_kabbandung[2, 1] / sum(freq_parsial_kabbandung[2, ])
RDKabBandung <- p1 - p2

#Z = Kota Cimahi
p1 <- freq_parsial_kotacimahi[1, 1] / sum(freq_parsial_kotacimahi[1, ])
p2 <- freq_parsial_kotacimahi[2, 1] / sum(freq_parsial_kotacimahi[2, ])
RDKotaCimahi <- p1 - p2

#Z = Kabupaten Bandung Barat
p1 <- freq_parsial_kabbandungbrt[1, 1] / sum(freq_parsial_kabbandungbrt[1, ])
p2 <- freq_parsial_kabbandungbrt[2, 1] / sum(freq_parsial_kabbandungbrt[2, ])
RDKabBandungBarat <- p1 - p2

list(Kota_Bandung = RDKotaBandung, Kab_Bandung = RDKabBandung, Kota_Cimahi = RDKotaCimahi, Kab_BandungBarat = RDKabBandungBarat)

## $Kota_Bandung
## [1] -0.01418648
## 
## $Kab_Bandung
## [1] -0.01487922
## 
## $Kota_Cimahi
## [1] -0.09247751
## 
## $Kab_BandungBarat
## [1] -0.01920472

Pada Kota Bandung, perempuan memiliki perempuan memiliki risiko 1.42% lebih tinggi dalam bekerja dibandingkan laki-laki. Pada Kabupaten Bandung, perempuan memiliki perempuan memiliki risiko 1.48% lebih tinggi dalam bekerja dibandingkan laki-laki. Pada Kota Cimahi, perempuan memiliki perempuan memiliki risiko 9.2% lebih tinggi dalam bekerja dibandingkan laki-laki. Pada Kabupaten Bandung Barat, perempuan memiliki perempuan memiliki risiko 3.2% lebih tinggi dalam bekerja dibandingkan laki-laki.

Risiko Relative (Relative Risk , RR)

\[ RR = \frac{P(Y|X1,Z)}{P(Y|X2,Z)} \]

Perhitungan menggunakan R :

#Z = Kota Bandung
p1 <- freq_parsial_kotabandung[1, 1] / sum(freq_parsial_kotabandung[1, ])
p2 <- freq_parsial_kotabandung[2, 1] / sum(freq_parsial_kotabandung[2, ])
RRKotaBandung <- p1 / p2

#Z = Kabupaten Bandung
p1 <- freq_parsial_kabbandung[1, 1] / sum(freq_parsial_kabbandung[1, ])
p2 <- freq_parsial_kabbandung[2, 1] / sum(freq_parsial_kabbandung[2, ])
RRKabBandung <- p1 / p2

#Z = Kota Cimahi
p1 <- freq_parsial_kotacimahi[1, 1] / sum(freq_parsial_kotacimahi[1, ])
p2 <- freq_parsial_kotacimahi[2, 1] / sum(freq_parsial_kotacimahi[2, ])
RRKotaCimahi <- p1 / p2

#Z = Kabupaten Bandung Barat
p1 <- freq_parsial_kabbandungbrt[1, 1] / sum(freq_parsial_kabbandungbrt[1, ])
p2 <- freq_parsial_kabbandungbrt[2, 1] / sum(freq_parsial_kabbandungbrt[2, ])
RRKabBandungBarat <- p1 / p2

list(Kota_Bandung = RRKotaBandung, Kab_Bandung = RRKabBandung, Kota_Cimahi = RRKotaCimahi, Kab_BandungBarat = RRKabBandungBarat)

## $Kota_Bandung
## [1] 0.9845889
## 
## $Kab_Bandung
## [1] 0.9842439
## 
## $Kota_Cimahi
## [1] 0.90282
## 
## $Kab_BandungBarat
## [1] 0.9793863

Pada Kota Bandung, laki-laki memliki risiko 0.9845 kali lebih kecil mendapatkan pekerjaan dibandingkan perempuan. Pada Kabupaten Bandung, laki-laki memliki risiko 0.9842 kali lebih kecil mendapatkan pekerjaan dibandingkan perempuan. Pada Kota Cimahi, laki-laki memliki risiko 0.9028 kali lebih kecil mendapatkan pekerjaan dibandingkan perempuan. Pada Kabupaten Bandung Barat, laki-laki memliki risiko 0.979 kali lebih kecil mendapatkan pekerjaan dibandingkan perempuan.

Odds Ratio (OR)

\[ OR = \frac{P(Y|X1,Z)/(1-P(Y|X1,Z))}{P(Y|X2,Z)/(1-P(Y|X2,Z))} \]

Perhitungan menggunakan R :

#Z = Kota Bandung
o1 <- freq_parsial_kotabandung[1, 1] * freq_parsial_kotabandung[2,2]
o2 <- freq_parsial_kotabandung[1, 2] * freq_parsial_kotabandung[2,1]
ORKotaBandung <- o1 / o2

#Z = Kabupaten Bandung
o1 <- freq_parsial_kabbandung[1, 1] * freq_parsial_kabbandung[2,2]
o2 <- freq_parsial_kabbandung[1, 2] * freq_parsial_kabbandung[2,1]
ORKabBandung <- o1 / o2

#Z = Kota Cimahi
o1 <- freq_parsial_kotacimahi[1, 1] * freq_parsial_kotacimahi[2,2]
o2 <- freq_parsial_kotacimahi[1, 2] * freq_parsial_kotacimahi[2,1]
ORKotaCimahi <- o1 / o2

#Z = Kabupaten Bandung Barat
o1 <- freq_parsial_kabbandungbrt[1, 1] * freq_parsial_kabbandungbrt[2,2]
o2 <- freq_parsial_kabbandungbrt[1, 2] * freq_parsial_kabbandungbrt[2,1]
ORKabBandungBarat <- o1 / o2

list(Kota_Bandung = ORKotaBandung, Kab_Bandung = ORKabBandung, Kota_Cimahi = ORKotaCimahi, Kab_BandungBarat = ORKabBandungBarat)

## $Kota_Bandung
## [1] 0.8354417
## 
## $Kab_Bandung
## [1] 0.7766072
## 
## $Kota_Cimahi
## [1] 0.3101283
## 
## $Kab_BandungBarat
## [1] 0.7645645

Pada Kota Bandung, Laki-Laki memiliki 0.83 lebih kecil odds mendapatkan pekerjaan dibandingkan Perempuan. Pada Kabupaten Bandung, Laki-Laki memiliki 0.77 lebih kecil odds mendapatkan pekerjaan dibandingkan Perempuan. Pada Kota Cimahi, Laki-Laki memiliki 0.31 lebih kecil odds mendapatkan pekerjaan dibandingkan Perempuan. Pada Kabupaten Bandung Barat, Laki-Laki memiliki 0.76 lebih kecil odds mendapatkan pekerjaan dibandingkan Perempuan.

Conditional Independence

Kondisi ini terjadi ketika X dan Y tidak saling bergantung setelah dikontrol oleh Z.

\[ P(X,Y|Z) = P(X|Z)P(Y|Z) \]

Pada kasus ini, Menguji apakah Gender dan Status Pekerjaan independen secara kondisional terhadap Wilayah.

Perhitungan menggunakan R :

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

# Uji Chi-Square untuk masing-masing wilayah
 chisq_kotabandung <- chisq.test(data2[,,"Kota Bandung"])
 chisq_kabbandung <- chisq.test(data2[,,"Kabupaten Bandung"])
 chisq_kotacimahi <- chisq.test(data2[,,"Kota Cimahi"])
 chisq_kabbandungbarat <- chisq.test(data2[,,"Kabupaten Bandung Barat"])
 
 list(Kota_Bandung = chisq_kotabandung, Kab_Bandung = chisq_kabbandung, Kota_Cimahi = chisq_kotacimahi, Kab_BandungBarat = chisq_kabbandungbarat)

## $Kota_Bandung
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data2[, , "Kota Bandung"]
## X-squared = 775.92, df = 1, p-value < 2.2e-16
## 
## 
## $Kab_Bandung
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data2[, , "Kabupaten Bandung"]
## X-squared = 1574.1, df = 1, p-value < 2.2e-16
## 
## 
## $Kota_Cimahi
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data2[, , "Kota Cimahi"]
## X-squared = 6795.5, df = 1, p-value < 2.2e-16
## 
## 
## $Kab_BandungBarat
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data2[, , "Kabupaten Bandung Barat"]
## X-squared = 1041.5, df = 1, p-value < 2.2e-16

Dari ke-empat daerah ini didapatkan bahwa p-value < \(\alpha\) (2.2e-16 < 0.05). Sehingga dapat disipulkan bahwa asih ada hubungan signifikan antara gender dan status pekerjaan meskipun wilayah dikendalikan

Marginal Y dan X

Menunjukkan distribusi gabungan Y dan X untuk setiap wilayah Z.

Z = Kota Bandung

Perhitungan Menggunakan R :

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

#Kota Bandung
a = data2["Laki-laki","Bekerja","Kota Bandung"]
b = data2["Perempuan","Bekerja","Kota Bandung"]
c = data2["Laki-laki","Tidak Bekerja","Kota Bandung"]
d = data2["Perempuan","Tidak Bekerja","Kota Bandung"]

p1 <- a / (a + b)
p2 <- c / (c + d)
RD_kotabandung <- p1- p2
RR_kotabandung <- p1 / p2
OR_kotabandung <- (a * d) / (b * c)
data_RD_RR_OR_kotabandung <- data.frame(RD = RD_kotabandung, RR = RR_kotabandung, OR = OR_kotabandung)
data_RD_RR_OR_kotabandung

Z = Kabupaten Bandung

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

#Kabupaten Bandung
a = data2["Laki-laki","Bekerja","Kabupaten Bandung"]
b = data2["Perempuan","Bekerja","Kabupaten Bandung"]
c = data2["Laki-laki","Tidak Bekerja","Kabupaten Bandung"]
d = data2["Perempuan","Tidak Bekerja","Kabupaten Bandung"]

p1 <- a / (a + b)
p2 <- c / (c + d)
RD_kabbandung <- p1- p2
RR_kabbandung <- p1 / p2
OR_kabbandung <- (a * d) / (b * c)
data_RD_RR_OR_kabbandung <- data.frame(RD = RD_kabbandung, RR = RR_kabbandung, OR = OR_kabbandung)
data_RD_RR_OR_kabbandung

Z = Kota Cimahi

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

#Kota Cimahi
a = data2["Laki-laki","Bekerja","Kota Cimahi"]
b = data2["Perempuan","Bekerja","Kota Cimahi"]
c = data2["Laki-laki","Tidak Bekerja","Kota Cimahi"]
d = data2["Perempuan","Tidak Bekerja","Kota Cimahi"]

p1 <- a / (a + b)
p2 <- c / (c + d)
RD_kotacimahi <- p1- p2
RR_kotacimahi <- p1 / p2
OR_kotacimahi <- (a * d) / (b * c)
data_RD_RR_OR_kotacimahi <- data.frame(RD = RD_kotacimahi, RR = RR_kotacimahi, OR = OR_kotacimahi)
data_RD_RR_OR_kotacimahi

Z = Kabupaten Bandung Barat

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

#Kabupaten Bandung Barat
a = data2["Laki-laki","Bekerja","Kabupaten Bandung Barat"]
b = data2["Perempuan","Bekerja","Kabupaten Bandung Barat"]
c = data2["Laki-laki","Tidak Bekerja","Kabupaten Bandung Barat"]
d = data2["Perempuan","Tidak Bekerja","Kabupaten Bandung Barat"]

p1 <- a / (a + b)
p2 <- c / (c + d)
RD_kabbandungbrt <- p1- p2
RR_kabbandungbrt <- p1 / p2
OR_kabbandungbrt <- (a * d) / (b * c)
data_RD_RR_OR_kabbandungbrt <- data.frame(RD = RD_kabbandungbrt, RR = RR_kabbandungbrt, OR = OR_kabbandungbrt)
data_RD_RR_OR_kabbandungbrt

Inferensi Tabel Kontingensi 3 Arah

Tujuan dari inferensi dalam tabel kontingensi tiga arah (2×2×K) adalah untuk mengevaluasi hubungan antara dua variabel (misalnya X dan Y) dengan mempertimbangkan pengaruh variabel ketiga (Z). Ini sangat penting dalam mencegah kesalahan interpretasi akibat variabel perancu.

Independensi Bersyarat

Dua variabel X dan Y dikatakan independen bersyarat terhadap Z jika dalam setiap strata Z, hubungan antara X dan Y tidak signifikan.

\(OR(X,Y∣Z)=1\)

Tidak ada pengaruh langsung antara X dan Y setelah mengontrol Z

Pengujian Statistik

Uji CMH digunakan untuk menguji apakah X dan Y tetap berhubungan setelah efek Z dikontrol.

Hipotesis nol (H₀): X dan Y independen bersyarat terhadap Z.
Hipotesis alternatif (H₁): Ada asosiasi antara X dan Y setelah mengontrol Z.

Statistik Uji :

\[ CMH = \frac{[\sum (a_k-E(a_k))]^2}{\sum Var(a_k)} \]

dengan,

\(a_k\) = frekeunsi sel atas kiri di strata ke-k

\(E(a_k)\) = Nilai harapan dari \(a_k\)

\(Var(a_k)\) = varians dari \(a_k\)

Perhitungan Menggunakan R :

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

ujicmh = mantelhaen.test(data2, correct = FALSE)
ujicmh

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  data2
## Mantel-Haenszel X-squared = 6418.4, df = 1, p-value < 2.2e-16
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.7311315 0.7421972
## sample estimates:
## common odds ratio 
##         0.7366436

Didapatkan bahwa p-value < \(\alpha\) (2.2e-16 < 0.05) sehingga dapat disimpulkan bahwa terdapat hubungan signifikan antara Jenis Kelamin dan Status Pekerjaan setelah mengontrol Wilayah.

Odds Ratio bersama

Jika dalam setiap tabel parsial nilai odds ratio (OR) memiliki arah yang sama dan relatif stabil, kita bisa menghitung odds ratio gabungan atau common odds ratio menggunakan metode Mantel-Haenszel:

\[ \theta MH=\frac{\sum(n11k\times n22k/n..k)}{\sum(n12k\times n21k/n..k)} \]

Perhitungan menggunakan R :

library(epitools)
# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

ujimh = mantelhaen.test(data2, correct = FALSE)
ujimh

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  data2
## Mantel-Haenszel X-squared = 6418.4, df = 1, p-value < 2.2e-16
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.7311315 0.7421972
## sample estimates:
## common odds ratio 
##         0.7366436

Didapatkan bahwa p-value < \(\alpha\) (2.2e-16 < 0.05) sehingga dapat disimpulkan bahwa terdapat hubungan signifikan antara Jenis Kelamin dan Status Pekerjaan setelah mengontrol Wilayah.

Uji Homogenitas OR dengan Breslow-Day

Dilakukan untuk Menguji apakah odds ratio antar strata konsisten (homogen).

Hipotesis

H₀: OR sama untuk semua strata
H₁: Setidaknya ada satu strata dengan OR berbeda.

Statistik Uji :

\[ X^2BD = \sum \frac{(aj-\hat aj)^2}{Var(aj|OR)} \]

dengan,

\(j\) : indeks untuk strata

\(a_h\) = frekuensi observasi sel a pada strata ke-h

\(\hat a_h\) = frekuensi harapan berdasarkan OR gabungan

\(Var(\hat a_h)\) = varians dari frekuensi harapan

Perhitungan Menggunakan 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

library(DescTools)

## Warning: package 'DescTools' was built under R version 4.3.3

# Membuat tabel kontingensi 3 arah
data2 = array(
  c(742536, 459949, 76725, 39705,
    1116497,639831, 84723, 37706,
    166488, 115910, 27298, 5894,
    575757, 294512, 55248, 21607),
  dim = c(2, 2, 4),
  dimnames = list(
    Jenis_Kelamin = c("Laki-laki", "Perempuan"),
    Status_Pekerjaan = c("Bekerja", "Tidak Bekerja"),Wilayah = c("Kota Bandung","Kabupaten Bandung","Kota Cimahi","Kabupaten Bandung Barat")
  )
)

bd_test = BreslowDayTest(data2)
bd_test

## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  data2
## X-squared = 4034.1, df = 3, p-value < 2.2e-16

Didapatkan bahwa p-value < \(\alpha\) (2.2e-16 < 0.05) sehingga rasio odds tidak homogen di semua wilayah.

GLM

GLM adalah perluasan dari model regresi linear klasik yang memungkinkan analisis terhadap data non-normal, seperti data biner dan cacah. GLM terdiri dari tiga komponen utama:

Distribusi dari exponential family
Fungsi link: menghubungkan rata-rata respons dengan kombinasi linear prediktor.
Fungsi linear prediktor:

\[\eta=X\beta\]

Exponential Family

Distribusi termasuk dalam exponential family jika dapat ditulis dalam bentuk:

\[f(y; \theta, \phi) = \exp\left\{\frac{y\theta - b(\theta)}{\phi} + c(y, \phi)\right\}\]

Contoh distribusi dalam keluarga ini:

Normal
Binomial
Poisson
Gamma

Keluarga ini memungkinkan estimasi dan inferensi statistik melalui pendekatan GLM.

Studi Kasus

Dari suatu pabrik wine / anggur di Italia, pemilik pabrik ingin menentukan apakakah suatu anggur termasuk ke dalam tipe 1, tipe 2, atau tipe 3 berdasarkan kadar alkohol dalam anggur (Alcohol). Peneliti ingin melihat apakah suatu wine termasuk tipe 1 atau buka Oleh karena itu, terdapat 2 kemungkinan peluang. Yaitu 0 untuk tipe 1 dan 1 untuk bukan tipe 1. Buatlah dalam Model Regresi Logistik dan Poisson untuk menentukan apakah wine termasuk tipe 1 atau bukan ?

Sumber Data : Aeberhard, S., and D. Coomans. “Wine.” UCI Machine Learning Repository, https://archive.ics.uci.edu/dataset/109/wine.

Implementasi dalam R :

#Input Data
getwd()

## [1] "D:/SEMESTER 4/ADK/E-BOOK UTS"

setwd("D:/SEMESTER 4/ADK/E-BOOK UTS")
data3 = read.csv("dataEBOOK2.csv", header = T, sep = ",")
data3

attach(data3)

Model Regresi Logistik

Model Regresi Logistik digunakan untuk memodelkan data biner (0 atau 1)

Keunggulan Utama Regresi Logistik:

Menghasilkan probabilitas (0–1), bukan hanya klasifikasi.
Tidak mengharuskan hubungan linear antara X dan Y, seperti pada regresi linear biasa.
Cocok untuk variabel prediktor campuran: numerik, kategorik, atau gabungan keduanya.
Interpretasi parameter mudah melalui odds ratio.
Mudah diperluas ke regresi logistik multikategori atau regresi logistik ordinal.

Fungsi Link (logit) :

\[ g(\mu) = log(\frac{\mu}{1-\mu}) \]

Model Regresi Logistik :

\[ log(\frac{\mu}{1-\mu})=X\beta \]

Fungsi Inverse Link :

\[ \mu = \frac{exp(X\beta)}{1+exp(X\beta)} \]

Estimasi Parameter

Metode estimasi parameter pada GLM umumnya menggunakan Maximum Likelihood Estimation (MLE). Log-likelihood fungsi untuk regresi logistik:

\[ l(\beta)=\sum [yi\ log(\pi i)+(1-yi)log(1-\pi i)] \]

dengan,

\[ \pi i=\frac{exp(xi^T\beta)}{1+exp(xi^T\beta)} \]

Estimasi dilakukan melalui iterasi Newton-Rapshon dengan perhitungan melalui R :

#Input Data
getwd()

## [1] "D:/SEMESTER 4/ADK/E-BOOK UTS"

setwd("D:/SEMESTER 4/ADK/E-BOOK UTS")
data3 = read.csv("dataEBOOK2.csv", header = T, sep = ",")
data3

modellog = glm(Class.Alcohol..Y. ~ Alcohol..X1., data = data3, family = binomial )
summary(modellog)

## 
## Call:
## glm(formula = Class.Alcohol..Y. ~ Alcohol..X1., family = binomial, 
##     data = data3)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -31.6592     6.7901  -4.663 3.12e-06 ***
## Alcohol..X1.   2.3685     0.5049   4.691 2.72e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 147.20  on 106  degrees of freedom
## Residual deviance: 115.03  on 105  degrees of freedom
## AIC: 119.03
## 
## Number of Fisher Scoring iterations: 4

Setelah dilakukan perhitungan didapatkan bahwa Model Regresi Logistik nya adalah :

\[ Y = -31.659 + 2.368 \times X1 \]Dengan kata lain, setiap kenaikan 1 persen alkohol akan meningkatkan 2.368 nilai Y dengan nilai awal -31.659

Selain itu p-value < \(\alpha\) (2.72e-06 < 0.05) didapatkan bahwa Alkohol memengaruhi Tipe dari Wine tersebut.

Visualisasi menggunakan plot :

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.3

 data3$pred <- predict(modellog, type = "response")
 ggplot(data3, aes(x =Alcohol..X1., y = Class.Alcohol..Y.)) +
 geom_point(alpha = 0.5, color = "gray40") +
 geom_line(aes(y = pred), color = "blue", linewidth = 1.5) +
 labs(title = "Kurva Logit pada Regresi Logistik",
 x="X (Prediktor)",
 y="Probabilitas / Respons") +
 theme_minimal()

Dari Kurva ini dapat terlihat bahwa Bentuk sigmoid ini menunjukkan hubungan antara prediktor dan probabilitas hasil bernilai 1.

Model Regresi Poisson

Model Regresi Poisson digunakan ketika variabel respon merupakan data cacah (count), seperti jumlah kecelakaan, panggilan telepon, kasus penyakit, dll.

\[ P(Y=y)=\frac{e^{-\lambda}\lambda^y}{y!} \]

Hal ini dapat ditulis dalam bentuk exponential family :

\[ f(y;\theta)=exp[y\ log(\lambda)-\lambda-log(y!)] \]

Fungsi link :

\[ g(\mu)=log(\mu) \]

Model Regresi :

\[ log(\mu i)=xi^T\beta \]

Fungsi Inverse Link :

\[ \mu i=exp(xi^T\beta) \]

Estimasi Parameter

Metode estimasi parameter pada GLM umumnya menggunakan Maximum Likelihood Estimation (MLE):

\[ l(\beta) = \sum[yi\ log(\lambda i)-\lambda i -log(yi!)] \]

dengan,

\[ \lambda i = exp(xi^T\beta) \]

Perhitungan dilakukan menggunakan iterasi Newton-Raphson dengan R sebagai berikut :

#Input Data
getwd()

## [1] "D:/SEMESTER 4/ADK/E-BOOK UTS"

setwd("D:/SEMESTER 4/ADK/E-BOOK UTS")
data3 = read.csv("dataEBOOK2.csv", header = T, sep = ",")
data3

modelpois = glm(Class.Alcohol..Y. ~ Alcohol..X1., data = data3, family = poisson )
summary(modelpois)

## 
## Call:
## glm(formula = Class.Alcohol..Y. ~ Alcohol..X1., family = poisson, 
##     data = data3)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -12.1849     3.3084  -3.683 0.000231 ***
## Alcohol..X1.   0.8513     0.2405   3.539 0.000401 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 70.244  on 106  degrees of freedom
## Residual deviance: 57.165  on 105  degrees of freedom
## AIC: 179.16
## 
## Number of Fisher Scoring iterations: 5

Setelah dilakukan perhitungan didapatkan bahwa Model Regresi Poisson nya adalah :

\[ Y = -12.1849 + 0.8512 \times X1 \] Dengan kata lain, setiap kenaikan 1 persen alkohol akan meningkatkan 0.8512 nilai Y dengan nilai awal -12.1849

Selain itu p-value < \(\alpha\) (0.000401 < 0.05) didapatkan bahwa Alkohol memengaruhi Tipe dari Wine tersebut.

Visualisasi menggunakan plot :

library(ggplot2)
 data3$pred <- predict(modelpois, type = "response")
 ggplot(data3, aes(x =Alcohol..X1., y = Class.Alcohol..Y.)) +
 geom_point(alpha = 0.5, color = "gray40") +
 geom_line(aes(y = pred), color = "blue", linewidth = 1.5) +
 labs(title = "Kurva Logit pada Regresi Logistik",
 x="X (Prediktor)",
 y="Probabilitas / Respons") +
 theme_minimal()

Dari kurva ini terlihat bahwa kurva naik eksponensial, yang artinya semakin besar X maka jumlah kejadian semakin besar.

Inferensi GLM

Dalam Generalized Linear Model (GLM), inferensi statistik membutuhkan pemahaman terhadap ekspektasi dan varians dari estimator model, terutama untuk mengembangkan alat-alat uji seperti Wald test, Likelihood Ratio test, dan interval kepercayaan.

Ekspektasi Estimator :

Estimator dikatakan tak bias jika nilai harapannya sama dengan parameter yang diperkirakan, yaitu:

\[ E[\hat{\beta}]= \beta \]

Varians Estimator

Varians digunakan untuk mengetahui presisi estimasi:

\[ Var(\hatβ)≈(X^⊤WX)^{−1} \]

dengan W adalah matriks bobot tergantung pada distribusi dan link function.

Varians dalam GLM tidak konstan, Tidak seperti regresi linear (OLS) yang mengasumsikan homoskedastisitas, Pada GLM berlaku :

\[ Var(Yi) =\phi V(\mu i) \]

Mencari Ekspektasi dan Varians dalam GLM

Ekspektasi :

Turunan pertama dari fungsi momen sebagai berikut :

\[ E(Y) = \int yf(y;\theta)dy = \mu \]

Ekspektasi Turunan Pertama :

\[ U(\theta) = \frac{\partial l}{\partial \theta}=y - b^\prime(\theta) \]

\[ E[U(\theta)]=E[y-b^\prime(\theta)]=\mu - b^\prime (\theta) = 0 \]

Maka didapatkan :

\[ \mu=b^\prime(\theta) \]

Varians :

Turunan kedua dari fungsi momen :

\[ \frac{\partial^2l}{\partial\theta^2}=-b^{\prime \prime}(\theta) \]

Sehingga didapatkan :

\[ Var(Y)=b^{\prime \prime}(\theta) = \phi V(\mu) \]

Metode Penaksiran Parameter

GLM menggunakan metode Maximum Likelihood Estimation (MLE). Dengan Turunan Pertama = 0 dan Turutan Kedua < 0. Namun karena bentuk GLM tidak eksplisit digunakan metode Optimisasi Newton-Raphson.

Iterasi :

\[ \beta^{(t+1)}=\beta^{(t)}-H^{-1}(\beta^{(t)})U(\beta^{(t)}) \]

Diagnostik Model GLM

Digunakan untuk mengevaluasi kecocokan model.

Statistik Devians :

Dilakukan untuk mengukur apakah ada model lain yang lebih baik ?

Devians :

\[ D = 2\sum[yi\ log(\frac{yi}{\hat\mu i})-(yi-\hat\mu i)] \]

Statistik Chi-Square Pearson :

Menguji apakah model lebih baik daripada tidak ada model sama sekali

Statistik :

\[ X^2=\sum \frac{(yi-\hat \mu i)^2}{\hat \mu i} \]

Detail Metode Estimasi dan Inferensi Regresi Logistik

Estimasi Parameter

Mencari parameter \(\beta\) yang memaksimlakan fungsi likelihood dengan Metode Newton-Rapshon

Model Regresi Logistik untuk probabilitas :

\[ \pi i= \frac{1}{1+exp(-xi^T\beta)} \]

Fungsi likelihood untuk n observasi :

\[ l(\beta)=\sum [yi\ log(\pi i)+(1-yi)log(1-\pi i)] \]

Langkah-Langkah :

Mencari Turunan Pertama
Mencari Turunan Kedua
Iterasi Newton-Raphson

Perhitungan Menggunakan R :

X <- cbind(1, Alcohol..X1.)  
beta <- matrix(c(0, 0), nrow = 2)


for (i in 1:10) {
  z <- X %*% beta
  p <- 1 / (1 + exp(-z))
  W <- diag(as.vector(p * (1 - p)))
  gradient <- t(X) %*% (Class.Alcohol..Y. - p)
  hessian <- -t(X) %*% W %*% X
  
  beta_new <- beta - solve(hessian) %*% gradient
  
  cat("Iterasi", i, ": beta =\n")
  print(beta_new)
  
  if (max(abs(beta_new - beta)) < 1e-6) break
  
  beta <- beta_new
}

## Iterasi 1 : beta =
##                    [,1]
##              -23.998491
## Alcohol..X1.   1.795606
## Iterasi 2 : beta =
##                    [,1]
##              -30.678582
## Alcohol..X1.   2.295182
## Iterasi 3 : beta =
##                    [,1]
##              -31.641259
## Alcohol..X1.   2.367145
## Iterasi 4 : beta =
##                    [,1]
##              -31.659154
## Alcohol..X1.   2.368482
## Iterasi 5 : beta =
##                    [,1]
##              -31.659160
## Alcohol..X1.   2.368482
## Iterasi 6 : beta =
##                    [,1]
##              -31.659160
## Alcohol..X1.   2.368482

Setelah dilakukan Iterasi menggunakan Newton-Raphson didapatkan hasil akhir estimasi sampai iterasi ke 6 dengan model seperti berikut :

\[ Y = -31.659160 + 2.368482\times X1 \]

Inferensi Parameter

Uji Wald

Uji Wald digunakan untuk menguji signifikansi parameter \(\beta j\) dalam model regresi logistik :

Statistik Uji :

\[ W =Z^2 = (\frac{\hat\beta j}{SE(\hat\beta j)})^2 ∼ \chi_i^2 \]

Perhitungan menggunakan R :

#Ambil Nilai Koef dan SE
beta_hat <- coef(modellog)["Alcohol..X1."]
se_beta <- summary(modellog)$coefficients["Alcohol..X1.", "Std. Error"]
 
#Hitung Z
Z = beta_hat / se_beta
Z

## Alcohol..X1. 
##     4.690908

Wald_stat = Z^2
Wald_stat

## Alcohol..X1. 
##     22.00461

p_value = 1 - pchisq(Wald_stat, df = 1)
p_value

## Alcohol..X1. 
## 2.719959e-06

Karena p-value < \(\alpha\) (2.719959e-06 < 0.05 ) didapatkan bahwa variabel Alcohol berpengaruh terhdap kelas alkohol dari suatu wine.

Uji Likelihood Ratio

Membandingkan model penuh dengan model tanpa prediktor.

Statistik Uji :

\[ G^2 = -2(log\ L_{null}-log\ L_{full})∼\chi^2 \]

Perhitungan menggunakan R :

# Model null
model_null <- glm(Class.Alcohol..Y.~ 1, data = data3, family = binomial)
# Likelihood ratio test
anova(model_null, model_null, test = "Chisq")

Evaluasi Kebaikan Model

Akaike Information Criterion (AIC)

Semakin kecil AIC, semakin baik model

Perhitungan menggunakan R :

AIC(modellog)

## [1] 119.0334

Bayesian Information Criterion (BIC)

Alternatif terhadap AIC, menghukum kompleksitas model

Perhitungan menggunakan R :

BIC(modellog)

## [1] 124.3791

Detail Metode Estimasi dan Inferensi Regresi Poisson

Estimasi Parameter

Estimasi ini juga dilakukan dengan memaksimalkan fungsi likelihood ratio dengan metode Itteratively Reweighted Least Squares (IRLS) seperti berikut :

Definisikan Model Regresi Poisson :

\[ log(\lambda i)=x_i^T\beta \]

sehingga, \(\lambda_i= exp (x_i^T\beta)\)
Mencari Log-Likelihood yang maksimal :

\[ l(\beta)=\sum[y_i\ log(\lambda_i)-\lambda_i-log(y_i!)] \]
Formulasi iteratif :

\[ \beta^{(t+1)}=(X^TW^{(t)}X)^{-1}X^TW^{(t)}z^{(t)} \]

dengan, \(W = diag(\lambda_i)\) , \(z=\eta + \frac{y-\lambda}{\lambda}\) , dan \(\eta_i=log(\lambda_i)=x_i^T\beta\)

Perhitungan Menggunakan R :

# Inisialisasi
 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 + (Class.Alcohol..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- 6

 beta

##                     [,1]
##              -12.1849446
## Alcohol..X1.   0.8512705

Setelah dilakukan Iterasi menggunakan IRLS didapatkan hasil akhir estimasi sampai iterasi ke 6 dengan model seperti berikut :

\[ Y = -12.1849446 + 0.8512705 \times X1 \]

Inferensi Parameter

Uji Wald

Uji Wald digunakan untuk menguji signifikansi parameter \(\beta j\) dalam model regresi logistik :

Statistik Uji :

\[ W =Z^2 = (\frac{\hat\beta j}{SE(\hat\beta j)})^2 ∼ \chi_i^2 \]

Perhitungan menggunakan R :

# Koefisien dan standar error
 coef_val <- coef(modelpois)[2]
 se_val <- summary(modelpois)$coefficients[2, 2]
 wald_z <- coef_val / se_val
 wald_chisq <- wald_z^2
 p_value <- 1- pchisq(wald_chisq, df = 1)
 cat("Z:", wald_z, "\nChi-Square:", wald_chisq, "\np-value:", p_value)

## Z: 3.539296 
## Chi-Square: 12.52662 
## p-value: 0.0004011954

Karena p-value < \(\alpha\) (0.0004011954< 0.05 ) didapatkan bahwa variabel Alcohol berpengaruh terhdap kelas alkohol dari suatu wine.

Uji Likelihood Ratio

Membandingkan model penuh dengan model tanpa prediktor.

Statistik Uji :

\[ G^2 = -2(log\ L_{null}-log\ L_{full})∼\chi^2 \]

Perhitungan menggunakan R :

# Model null 
model_null <- glm(Class.Alcohol..Y.~ 1, data = data3, family = poisson) 
# Likelihood ratio test 
anova(model_null, model_null, test = "Chisq")

Evaluasi Kebaikan Model

Akaike Information Criterion (AIC)

Semakin kecil AIC, semakin baik model

Perhitungan menggunakan R :

AIC(modelpois)

## [1] 179.1649

Bayesian Information Criterion (BIC)

Alternatif terhadap AIC, menghukum kompleksitas model. Perhitungan menggunakan R :

BIC(modelpois)

## [1] 184.5106

Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio

Regresi logistik digunakan untuk memodelkan hubungan antara variabel respon biner (misalnya: “Sakit” vs “Sehat”) dengan satu atau lebih variabel prediktor (X).

Dalam penerapannya, jenis skala pengukuran prediktor sangat memengaruhi cara kita memasukkannya ke dalam model.

Skala Nominal

Ciri: Tidak punya urutan logis; hanya membedakan kategori.
Contoh: Jenis Kelamin (Laki-laki, Perempuan), Warna, Wilayah.

Cara Masuk ke Model:

Harus diubah ke variabel dummy.
Misal:

-    Jenis Kelamin (Male = 0, Female = 1)

-    Wilayah: Barat, Timur, Utara → buat 2 dummy:

    -    Timur = 1 jika Timur, 0 lainnya

    -    Utara = 1 jika Utara, 0 lainnya (Barat sebagai referensi)

Skala Ordinal

Ciri: Ada urutan antar kategori, tapi jarak antar level belum tentu sama.
Contoh: Tingkat Pendidikan, Skor Kepuasan (1 = tidak puas, …, 5 = sangat puas).
Cara Masuk ke Model:

-   **Opsional**, tergantung konteks:

    1.   **Sebagai nominal** → Buat dummy untuk tiap level (seperti nominal)

    2.   📈 **Sebagai numerik** → Kodekan sebagai angka berurutan (misalnya: SMA = 1, S1 = 2, S2 = 3), **jika mengasumsikan efek linear atau hampir linear antar level**

Skala Rasio

Ciri: Numerik, punya jarak dan nol absolut (rasio bermakna).
Contoh: Pendapatan, Umur, Waktu.
Cara Masuk ke Model:

-   Bisa **langsung digunakan** sebagai prediktor numerik dalam model logistik.

-   Bisa juga **dikonversi ke kategorikal** jika diperlukan, misalnya pendapatan → kelompok (rendah, sedang, tinggi)

Simulasi Data

Suatu Penelitian ingin meneliti seberapa besar pengaruh Jenis Kelamin, Usia, dan Tingkat Pendidikan terhadap persepsi lamanya pelayanan Rumah Sakit Hasan Sadikin. Survei ini dilakukan pada tanggal 3 - 4 Juni 2025 dan dilakukan pada Ruang Tunggu IGD, ICU, Loket Pengambilan Hasil Lab, dan Radiologi.

X1 : Jenis Kelamin (Laki/Perempuan)

X2 : Tingkat Pendidikan (SD, SMP, SMA, Tinggi)

X3 : Usia

Y = Kecepatan Pelayanan Rumah Sakit (0 = Lama / Gagal, 1 = Cepat / Sukses)

library(tibble)

## Warning: package 'tibble' 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:vcdExtra':
## 
##     summarise

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# Simpan langsung hasil data yang dimutasi ke sim_data
sim_data <- tibble(
  JenisKelamin = c("Perempuan","Laki","Laki","Perempuan","Perempuan","Laki","Laki","Laki","Perempuan",
                   "Laki","Laki","Perempuan","Laki","Laki","Laki","Perempuan","Laki","Perempuan",
                   "Laki","Perempuan","Laki","Perempuan","Perempuan","Perempuan","Laki","Perempuan",
                   "Perempuan","Perempuan","Perempuan","Laki","Laki","Perempuan","Perempuan","Laki",
                   "Laki","Laki","Laki","Laki","Laki","Perempuan","Laki","Perempuan","Laki","Perempuan",
                   "Laki","Laki","Laki","Laki","Laki","Perempuan","Perempuan","Laki","Perempuan",
                   "Perempuan","Perempuan","Perempuan","Laki","Laki","Laki"),
  Pendidikan = c("SMP","SMA","Tinggi","SMP","Tinggi","SMA","SMA","SMA","Tinggi",
                 "SMA","SMA","SMA","SMA","SD","Tinggi","Tinggi","SMA","Tinggi",
                 "Tinggi","SMP","SMA","SMA","Tinggi","SMP","SMP","Tinggi","SD",
                 "SMA","SMP","SMA","SMA","Tinggi","SMA","SMA","Tinggi","Tinggi",
                 "Tinggi","SMA","SMA","SMA","SMA","SMA","SD","SMA","Tinggi",
                 "SMA","SMA","SMA","SMA","SMP","SMA","SMP","SMP","SMP",
                 "SMA","SD","Tinggi","SMA","SMA"),
  Umur = c(35,54,32,33,26,23,30,39,27,
           54,42,21,57,57,42,53,18,38,
           51,32,28,19,36,48,29,50,50,
           39,45,32,57,51,45,21,53,39,
           31,52,46,29,49,49,48,49,20,
           37,64,24,23,25,32,28,47,35,
           34,29,47,39,40),
  ResponAsli = c(4,4,3.5,4,4,3,4,5,4,
                 3,5,5,3,5,2,4,4,3,
                 1,4,3,4,3,5,5,3,5,
                 1,4,3,4,4,4,3,4,3,
                 4,4,4,3,3,1,3,3,4,
                 2,3,3,3,4,4,3.5,5,5,
                 4,5,4,3,3)
) %>%
  mutate(
    ResponCepat = case_when(
      ResponAsli >= 4 ~ 1,
      ResponAsli < 4 ~ 0,
      TRUE ~ NA_real_
    )
  ) %>%
  select(JenisKelamin, Pendidikan, Umur, ResponCepat)

# Lihat hasil awal
head(sim_data)

Dataset berisi 59 observasi dengan variabel Jenis Kelamin, Pendidikan, Umur, dan Status Keberhasilan

Eksplorasi Data

sim_data %>%
dplyr::group_by(ResponCepat) %>%
dplyr::summarise(Jumlah = dplyr::n(),
Rata2_Usia = mean(Umur)
)

Sebanyak 26 pengunjung menganggap pelayanan di RSHS masih lama dengan rata-rata usia responden yang menjawab adalah 39.269 tahun. Sedankgan 33 pengunjung menganggap pelayanan di RSHS sudah cepat dengan rata-rata usia responden yang menjawab adalah 38.242 tahun.

Perlakuan Variabel Ordinal

1. Treat sebagai Nominal (Dummy)

library(dplyr)
sim_data_nominal <- sim_data %>%
mutate(
education = factor(Pendidikan, levels = c("SD", "SMP", "SMA", "Tinggi"))
)
model_nominal <- glm(ResponCepat ~ JenisKelamin + education + Umur, data = sim_data_nominal, family = binomial)
summary(model_nominal)

## 
## Call:
## glm(formula = ResponCepat ~ JenisKelamin + education + Umur, 
##     family = binomial, data = sim_data_nominal)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)
## (Intercept)            0.840498   1.695541   0.496    0.620
## JenisKelaminPerempuan  0.942632   0.606682   1.554    0.120
## educationSMP           0.778981   1.627620   0.479    0.632
## educationSMA          -1.272382   1.262912  -1.007    0.314
## educationTinggi       -1.006892   1.309399  -0.769    0.442
## Umur                  -0.003305   0.024671  -0.134    0.893
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 80.959  on 58  degrees of freedom
## Residual deviance: 70.242  on 53  degrees of freedom
## AIC: 82.242
## 
## Number of Fisher Scoring iterations: 4

Keterangan: Model menggunakan education sebagai variabel dummy, baseline “HighSchool”. Setiap koefsien membandingkan kategori terhadap baseline.

Interpretasi Koefisien

Intercept (+1.04475)

Ini adalah log-odds dasar untuk individu berjenis kelamin Laki-Laki, berpendidikan SD, dan memiliki Usia = 0.
Tidak signifkan (p = 0.542), artinya baseline ini tidak berbeda secara signifkan dari probabilitas 50%.

Jenis Kelamin Perempuan(+0.889)

Individu Perempuan memiliki log-odds sukses yang lebih tinggi sebesar 0.889 dibandingkan Laki-Laki
Tidak signifkan (p = 0.144), artinya baseline ini tidak berbeda secara signifkan dari probabilitas 50%.

Education SMP (+0.7445)

Individu dengan Pendidikan SMP memiliki log-odds sukses yang lebih tinggi sebesar 0.7445 dibandingkan dengan Pendidikan SD.
Tidak signifkan (p = 0.647), artinya baseline ini tidak berbeda secara signifkan dari probabilitas 50%.

Education SMA (-1.2559)

Individu dengan Pendidikan SMA memiliki log-odds sukses yang lebih rendah sebesar -1,2559 dibandingkan dengan Pendidikan SD.
Tidak signifkan (p = 0.320), artinya baseline ini tidak berbeda secara signifkan dari probabilitas 50%.

Education Tinggi (-1.034)

Individu dengan Pendidikan Tinggi memiliki log-odds sukses yang lebih rendah sebesar -1,034 dibandingkan dengan Pendidikan SD.
Tidak signifkan (p = 0.430), artinya baseline ini tidak berbeda secara signifkan dari probabilitas 50%.

Umur (-0.007)

Setiap kenaikin 1 unit umur, mengurangi log-odds sebesar -0.007.
Tidak signifkan (p = 0.775), artinya umur tidak berhubungan positif dengan peluang cepat pelayanan.

Interpretasi Goodness-of Fit

Null deviance (80.959): Deviance model tanpa prediktor.
Residual deviance (70.242): Deviance model dengan prediktor. – Penurunan dari null deviance ke residual deviance menunjukkan model membawa informasi yang cukup baik.
AIC (82.242): – Digunakan untuk membandingkan model: semakin kecil AIC, semakin baik model dalam menyeimbangkan goodness-of-ft dan kompleksitas.

Signifkansi Model

Tidak ada Variabel yang Signifikan

Kesimpulan Praktis

Tidak ada Variabel yang bisa menjadi prediktor kuat seseorang akan menganggap pelayanan di RSHS cepat.

2. Treat sebagai Rasio (Numeric Rank)

sim_data_numeric <- sim_data %>%
mutate(
education_numeric = case_when(
Pendidikan == "SD" ~ 1,
Pendidikan == "SMP" ~ 2,
Pendidikan == "SMA" ~ 3,
Pendidikan == "Tinggi" ~ 4
)
)
model_numeric <- glm(ResponCepat ~ JenisKelamin + education_numeric + Umur, data = sim_data_numeric, family = binomial)
summary(model_numeric)

## 
## Call:
## glm(formula = ResponCepat ~ JenisKelamin + education_numeric + 
##     Umur, family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)  
## (Intercept)            1.411580   1.503316   0.939    0.348  
## JenisKelaminPerempuan  1.223990   0.577316   2.120    0.034 *
## education_numeric     -0.502749   0.359077  -1.400    0.161  
## Umur                  -0.004674   0.023781  -0.197    0.844  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 80.959  on 58  degrees of freedom
## Residual deviance: 73.179  on 55  degrees of freedom
## AIC: 81.179
## 
## Number of Fisher Scoring iterations: 4

Interpretasi : Model menggunakan education sebagai angka bertingkat 1-4, Koefisien education_numeric mengukur efek kenaikan satu tingkat pendidikan terhadap peluang sukses.

Interpretasi Koefisien

Intercept (+1.41158)

Ini adalah log-odds dasar untuk individu berjenis kelamin Laki-Laki, berpendidikan SD (education_numeric = 1) , dan memiliki Usia = 0.
Tidak signifkan (p = 0.348), artinya baseline ini tidak berbeda secara signifkan dari probabilitas 50%.

Jenis Kelamin Perempuan(+1.2239)

Individu Perempuan memiliki log-odds sukses yang lebih tinggi sebesar 1.2239 dibandingkan Laki-Laki
Signifkan (p = 0.034), artinya baseline ini berbeda secara signifkan dari probabilitas 50%.

Education Numeric (-0.5027)

Setiap kenaikan satu tingkat pendidikan (SD-> SMP -> SMA-> Tinggi) mengurangi log-odds sukses sebesar -0.5027.
Tidak signifkan (p = 0.161), artinya baseline ini tidak berbeda secara signifkan dari probabilitas 50%.

Umur (-0.0046)

Setiap kenaikin 1 unit umur, mengurangi log-odds sebesar -0.0046.
Tidak signifkan (p = 0.844), artinya umur tidak berhubungan positif dengan peluang cepat pelayanan.

Interpretasi Goodness-of Fit

Null deviance (80.959): Deviance model tanpa prediktor.
Residual deviance (73.179): Deviance model dengan prediktor. – Penurunan dari null deviance ke residual deviance menunjukkan model membawa informasi yang cukup baik.
AIC (81.179): – Digunakan untuk membandingkan model: semakin kecil AIC, semakin baik model dalam menyeimbangkan goodness-of-ft dan kompleksitas.

Signifkansi Model

Variabel Education Numeric signifikan terhadap kecepatan pelayanyan di RSHS
Variabel lainnya tidak signifikan

Kesimpulan Praktis

Semakin tinggi edukasi seseorang semakin rendah pendapat orang tersebut mengenai kecepatan pelayanan di RSHS
Variabel lainnya tidak signifikan terhadap pendapat orang tersebut mengenai kecepatan pelayanan di RSHS

Perbandingan Model

Menggunakan AIC

list(
AIC_Nominal = AIC(model_nominal),
AIC_Numeric = AIC(model_numeric)
)

## $AIC_Nominal
## [1] 82.24189
## 
## $AIC_Numeric
## [1] 81.17875

Interpretasi : Model dengan nilai AIC lebih kecil lebih disukai. Membandingkan mana yang lebih baik:treat education sebagai nominal atau numeric. Dari AIC ini terlihat bahwa Model Numerik lebih baik dipakai

Menggunakan Goodness of Fit Model

nullmod <- glm(ResponCepat ~ 1, data = sim_data, 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.132376 (df=6)
## 
## $McFadden_R2_Numeric
## 'log Lik.' 0.09610004 (df=4)

Interpretasi : Dari Uji McFadden Rsquared terlihat bahwa Model Nominal memiliki nilai yang lebih tinggi dibandingkan model numeric. Sehingga Model Nominal merupakan model yang lebih baik dibandingkan model numeric.

Visualisasi Prediksi :

library(ggplot2)
sim_data_nominal <- sim_data_nominal %>% mutate(predicted = predict(model_nominal, type = "response"))
sim_data_numeric <- sim_data_numeric %>% mutate(predicted = predict(model_numeric, type = "response"))

# Plot untuk model nominal
sim_data_nominal %>%
  ggplot(aes(x = Umur, y = predicted, color = education)) +
  geom_point(alpha = 0.7, size = 3) +
  scale_color_brewer(palette = "Set2", name = "Pendidikan") +
  labs(
    title = "Prediksi Probabilitas Respon Cepat Nominal (Ordinal sbg Nominal)",
    x = "Umur",
    y = "Probabilitas Prediksi"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

sim_data_numeric %>%
  ggplot(aes(x = Umur, y = predicted, color = as.factor(education_numeric))) +
  geom_point(alpha = 0.7, size = 3) +
  scale_color_brewer(palette = "Set2", name = "Pendidikan") +
  labs(
    title = "Prediksi Probabilitas Respon Cepat Numerik (Ordinal sebagai Numeric)",
    x = "Umur",
    y = "Probabilitas Prediksi"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

Ringkasan Koefsien Model

library(broom)

## Warning: package 'broom' was built under R version 4.3.3

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

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

# Ringkasan Model Nominal
tidy(model_nominal) %>%
  kable(format = "simple", booktabs = TRUE, caption = "Ringkasan Koefisien Model dengan Ordinal sebagai Nominal") %>%
  kable_styling(latex_options = c("hold_position", "striped"))

## Warning in kable_styling(., latex_options = c("hold_position", "striped")):
## Please specify format in kable. kableExtra can customize either HTML or LaTeX
## outputs. See https://haozhu233.github.io/kableExtra/ for details.

Ringkasan Koefisien Model dengan Ordinal sebagai Nominal
term	estimate	std.error	statistic	p.value
(Intercept)	0.8404983	1.6955412	0.4957109	0.6200984
JenisKelaminPerempuan	0.9426322	0.6066822	1.5537496	0.1202442
educationSMP	0.7789808	1.6276198	0.4786012	0.6322224
educationSMA	-1.2723816	1.2629119	-1.0074983	0.3136954
educationTinggi	-1.0068925	1.3093990	-0.7689730	0.4419093
Umur	-0.0033053	0.0246707	-0.1339755	0.8934220

# Ringkasan Model Numeric
tidy(model_numeric) %>%
  kable(format = "simple", booktabs = TRUE, caption = "Ringkasan Koefisien Model dengan Ordinal sebagai Numeric") %>%
  kable_styling(latex_options = c("hold_position", "striped"))

## Warning in kable_styling(., latex_options = c("hold_position", "striped")):
## Please specify format in kable. kableExtra can customize either HTML or LaTeX
## outputs. See https://haozhu233.github.io/kableExtra/ for details.

Ringkasan Koefisien Model dengan Ordinal sebagai Numeric
term	estimate	std.error	statistic	p.value
(Intercept)	1.4115801	1.5033156	0.9389779	0.3477421
JenisKelaminPerempuan	1.2239898	0.5773156	2.1201398	0.0339943
education_numeric	-0.5027492	0.3590769	-1.4001155	0.1614787
Umur	-0.0046737	0.0237806	-0.1965320	0.8441938

Pemilihan Model Regresi Logistik dan Evaluasi

Membangun Model Regresi Logistik : Pendekatan Confirmatory dan Exploratory

Dalam analisis regresi logistik, pemilihan model sangat penting untuk mendapatkan model yang akurat dalam memprediksi probabilitas kejadian suatu peristiwa (respon biner: 0 atau 1). Terdapat dua pendekatan utama:

1. Confirmatory Approach (Pendekatan Konfirmatori)

Definisi
Digunakan ketika peneliti sudah memiliki teori atau hipotesis yang jelas mengenai hubungan antara variabel prediktor dan respon.

Ciri-ciri:

Model dibangun berdasarkan teori atau hasil penelitian sebelumnya.
Fokus utama: menguji efek yang diasumsikan signifikan, bukan mencari model terbaik secara statistik.
Model lengkap (full model) biasanya dibangun terlebih dahulu.
Pengujian dilakukan dengan membandingkan model dengan dan tanpa suatu efek, menggunakan uji seperti:
- Likelihood Ratio Test (LRT)
- Wald Test
Tidak dilakukan seleksi variabel otomatis seperti AIC atau stepwise.
Contoh:
Jika teori menyatakan bahwa variabel x1 dan x2 memengaruhi keputusan membeli produk, maka model dibangun dengan x1 dan x2, lalu diuji apakah x2 signifikan terhadap respon.

2. Exploratory Approach (Pendekatan Eksploratori)

Definisi
Digunakan ketika belum ada teori kuat yang mendasari, atau ketika peneliti ingin mengeksplorasi hubungan potensial antar variabel prediktor.

Ciri-ciri:

Model dibangun bertahap, berdasarkan informasi dari data.
Seleksi variabel otomatis digunakan, berdasarkan kriteria statistik seperti:
- Akaike Information Criterion (AIC)
- Deviance
- Log-likelihood
Tiga metode seleksi umum:
- Forward Selection: mulai dari model kosong → tambah variabel satu per satu jika signifikan.
- Backward Elimination: mulai dari model penuh → buang variabel tidak signifikan.
- Stepwise Selection: gabungan forward dan backward.
Tujuan:

Menemukan model yang parsimonious: cukup sederhana, tetapi mampu menjelaskan data dengan baik.

Pemilihan antara pendekatan Confrmatory dan Exploratory bergantung pada tujuan penelitian. Jika ingin menguji hipotesis tertentu, gunakan pendekatan Confrmatory. Jika ingin menemukan model terbaik berdasarkan data, gunakan pendekatan Exploratory. Dalam praktiknya, kedua pendekatan ini sering digunakan secara komplementer: teori digunakan sebagai dasar, dan seleksi eksploratori dilakukan untuk menyempurnakan model

Simulasi Data

Prediksi Kepatuhan Pasien Mengonsumsi Obat

Latar belakang:
Sebuah rumah sakit ingin mengetahui faktor-faktor yang memengaruhi kepatuhan pasien dalam mengonsumsi obat secara teratur selama masa pengobatan.

📄 Deskripsi Variabel:

Variabel	Keterangan
`y` (respon)	Kepatuhan minum obat (1 = Patuh, 0 = Tidak patuh)
`x1`	Tingkat kecemasan pasien (skor dari tes psikologis, skala normal, numerik)
`x2`	Jenis kelamin pasien (0 = Laki-laki, 1 = Perempuan)
`x3`	Jumlah kunjungan ke dokter selama bulan terakhir (numerik)

library(knitr)
library(dplyr)
library(ggplot2)
library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(caret)

## Warning: package 'caret' was built under R version 4.3.3

## Loading required package: lattice

## 
## Attaching package: 'lattice'

## The following object is masked from 'package:gnm':
## 
##     barley

## 
## Attaching package: 'caret'

## The following objects are masked from 'package:DescTools':
## 
##     MAE, RMSE

library(pROC)

## Warning: package 'pROC' was built under R version 4.3.3

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

library(DescTools)
set.seed(123)

n <- 200
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n, mean = 1)
lin_pred <- -0.3 + 1.5 * x1 - 1.2 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
df <- data.frame(y = as.factor(y), x1, x2, x3)
head(df)

Pemilihan Model Model Full

model_full <- glm(y ~ x1 + x2 + x3, data = df, family = binomial)
summary(model_full)

## 
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.4821     0.3009  -1.602   0.1091    
## x1            1.4830     0.2395   6.191 5.96e-10 ***
## x2           -0.7757     0.3489  -2.223   0.0262 *  
## x3            0.7102     0.1814   3.914 9.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.54  on 199  degrees of freedom
## Residual deviance: 205.47  on 196  degrees of freedom
## AIC: 213.47
## 
## Number of Fisher Scoring iterations: 4

Metode Stepwise: Forward, Backward, dan Kedua Arah

null_model <- glm(y ~ 1, data = df, family = binomial)
step_forward <- step(null_model, direction = "forward", scope = formula(model_full), trace = FALSE)
step_backward <- step(model_full, direction = "backward", trace = FALSE)
step_both <- step(null_model, direction = "both", scope = formula(model_full), trace = FALSE)
AIC(model_full, step_forward, step_backward, step_both)

Evaluasi Model: ROC dan AUC

pred_prob <- predict(step_both, type = "response")
roc_obj <- roc(df$y, pred_prob)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

plot(roc_obj,main = "Kurva ROC", col="green")

auc(roc_obj)

## Area under the curve: 0.8236

Pseudo R-Squared

PseudoR2(step_both, which = c("CoxSnell", "Nagelkerke", "McFadden"))

##   CoxSnell Nagelkerke   McFadden 
##  0.2990559  0.3992214  0.2569821

Tabel Klasifkasi dan Evaluasi

pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
conf_matrix <- confusionMatrix(factor(pred_class), df$y, positive = "1")
conf_matrix

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 83 29
##          1 23 65
##                                           
##                Accuracy : 0.74            
##                  95% CI : (0.6734, 0.7993)
##     No Information Rate : 0.53            
##     P-Value [Acc > NIR] : 8.656e-10       
##                                           
##                   Kappa : 0.4762          
##                                           
##  Mcnemar's Test P-Value : 0.4881          
##                                           
##             Sensitivity : 0.6915          
##             Specificity : 0.7830          
##          Pos Pred Value : 0.7386          
##          Neg Pred Value : 0.7411          
##              Prevalence : 0.4700          
##          Detection Rate : 0.3250          
##    Detection Prevalence : 0.4400          
##       Balanced Accuracy : 0.7373          
##                                           
##        'Positive' Class : 1               
##

conf_matrix$byClass[c("Sensitivity", "Specificity")]

## Sensitivity Specificity 
##   0.6914894   0.7830189

Metode Perbandingan Model dalam Regresi Logistik”

Tujuan Dokumen ini menyajikan cara membandingkan model regresi logistik menggunakan ukuran:

- Deviance

-AIC (Akaike Information Criterion)

- Likelihood-Ratio serta menjelaskan prinsip Parsimony dalam pemilihan model

library(MASS)
library(broom)
library(DescTools)
set.seed(123)

Pembuatan Model

model1 <- glm(y ~ x1, data = df, family = binomial)
model2 <- glm(y ~ x1 + x2, data = df, family = binomial)
model3 <- glm(y ~ x1 + x2 + x3, data = df, family = binomial)

Perbandingan AIC dan Deviance

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

Likelihood-Ratio Test

anova(model1, model2, test = "LRT")

anova(model2, model3, test = "LRT")

Prinsip Parsimony

Model yang lebih kompleks cenderung memiliki nilai AIC dan deviance yang lebih kecil, karena lebih mampu menyesuaikan diri terhadap data. Namun, kompleksitas berlebihan dapat menyebabkan overfitting.

Maka digunakan prinsip parsimony (kesederhanaan):

Pilih model yang cukup baik menjelaskan data, tidak terlalu rumit, dan mudah diinterpretasikan.

Rumus dan Penjelasan Ukuran Evaluasi Model

1. Akaike Information Criterion (AIC)

Rumus:

AIC= = -2 log L + 2k

log L: log-likelihood model.
k: jumlah parameter dalam model.

Penjelasan:

AIC menggabungkan kecocokan model dan kompleksitas model.
Semakin kecil AIC, model dianggap lebih baik.
AIC memberi penalti untuk model dengan terlalu banyak parameter (kompleksitas tinggi).

Prinsip Parsimony dalam AIC:

Jika dua model memiliki AIC yang berbeda sedikit (misal < 2), pilih model yang lebih sederhana.

2. Deviance

Rumus:

D = −2[log⁡Lmodel−log⁡Lmodel saturasi]

Model saturasi adalah model yang cocok 100% dengan data (fit sempurna).

Penjelasan:

Deviance mengukur sejauh mana model berbeda dari model sempurna.
Deviance kecil → model lebih baik karena mendekati model saturasi.

3. Likelihood Ratio Test (LRT)

Rumus:

G^2 = -2 (Log L0 - Log L1)

Penjelasan:

LRT digunakan untuk menguji apakah model kompleks secara signifikan lebih baik.
Diuji menggunakan distribusi chi-square.
Jika:
G^2 besar dan p-value < 0.05
→ Model kompleks lebih unggul secara statistik.

Evaluasi Tabel Klasifkasi dan Akurasi Model

pred_prob <- predict(model3, type = "response")
pred_class <- factor(ifelse(pred_prob >= 0.5, 1, 0))
conf_matrix <- confusionMatrix(pred_class, df$y, positive = "1")
conf_matrix

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 83 29
##          1 23 65
##                                           
##                Accuracy : 0.74            
##                  95% CI : (0.6734, 0.7993)
##     No Information Rate : 0.53            
##     P-Value [Acc > NIR] : 8.656e-10       
##                                           
##                   Kappa : 0.4762          
##                                           
##  Mcnemar's Test P-Value : 0.4881          
##                                           
##             Sensitivity : 0.6915          
##             Specificity : 0.7830          
##          Pos Pred Value : 0.7386          
##          Neg Pred Value : 0.7411          
##              Prevalence : 0.4700          
##          Detection Rate : 0.3250          
##    Detection Prevalence : 0.4400          
##       Balanced Accuracy : 0.7373          
##                                           
##        'Positive' Class : 1               
##

Sensitivitas dan Spesifsitas

• Sensitivitas: Kemampuan model mendeteksi kelas positif secara benar (True Positive Rate) \(Sensitivity = \frac{TP}{TP+FN}\) • Spesifsitas: Kemampuan model mendeteksi kelas negatif secara benar (True Negative Rate) \(Specitifity = \frac{TN}{TN+FP}\)

conf_matrix$byClass[c("Sensitivity", "Specificity")]

## Sensitivity Specificity 
##   0.6914894   0.7830189

Kesimpulan • Deviance yang kecil menunjukkan kecocokan model yang lebih baik. • AIC yang rendah menunjukkan keseimbangan antara kecocokan dan kompleksitas. • Likelihood Ratio Test mengevaluasi apakah model kompleks secara signifkan lebih baik. • Tabel klasifkasi membantu menilai kinerja prediksi aktual vs prediksi model. • Prinsip Parsimony mengutamakan model sederhana jika performanya mirip.

Detail ROC - Penjelasan Kurva ROC (Receiver Operating Characteristic)

1. Definisi Kurva ROC

Sumbu Y (True Positive Rate / Sensitivitas):

\(Sensitivity= True Positive Rate = \frac{TP}{TP+FN}\)
Sumbu X (False Positive Rate / 1 - Spesifisitas):

\(1-Specificity= False Positve Rate = \frac{FP}{FP+TN}\)

Garis Diagonal: Menunjukkan performa acak (random guess).
Kurva Baik: Mendekati pojok kiri atas, menunjukkan model yang kuat.

2. Pergerakan Kurva ROC terhadap Threshold (Cut-off)

Menurunkan Threshold:

-   Sensitivitas meningkat

-    Spesifisitas menurun

Menaikkan Threshold:

-    Sensitivitas menurun

-    Spesifisitas meningkat

3. Karakteristik Kurva ROC Ideal

Naik secara vertikal hingga Sensitivity=1
Lanjut secara horizontal menuju 1 - Specificity=1
AUC (Area Under the Curve) mendekati 1

4. Interpretasi AUC (Area Under the Curve)

Nilai AUC	Interpretasi
0.5	Tidak lebih baik dari tebak acak
> 0.7	Model cukup baik
> 0.9	Model sangat baik

AUC juga dikenal sebagai Concordance Index
→ Probabilitas bahwa model memberi skor lebih tinggi pada kasus positif dibanding negatif.

5. Kegunaan Kurva ROC

Membandingkan beberapa model klasifikasi
Menentukan threshold optimal:
- Apakah False Positive atau False Negative lebih penting?
- Contoh: Dalam deteksi penyakit, False Negative sering lebih fatal.

6. Visualisasi dalam R

plot(roc_obj)

auc(roc_obj)

## Area under the curve: 0.8236

7. Simulasi Pemilihan Threshold Optimal

pred = predict(model3, type= "response")
thresholds <- seq(0.1, 0.9, by = 0.05)
results <- data.frame(Threshold = thresholds)
results$Sensitivity <- sapply(thresholds, function(t) {
pred_class <- ifelse(pred >= t, 1, 0)
cm <- table(Pred = pred_class, Obs = df$y)
TP <- cm["1", "1"]
FN <- cm["0", "1"]
TP / (TP + FN)
})
results$Specificity <- sapply(thresholds, function(t) {
pred_class <- ifelse(pred >= t, 1, 0)
cm <- table(Pred = pred_class, Obs = df$y)
TN <- cm["0", "0"]
FP <- cm["1", "0"]
TN / (TN + FP)
})
print(results)

##    Threshold Sensitivity Specificity
## 1       0.10   1.0000000   0.1603774
## 2       0.15   0.9893617   0.3113208
## 3       0.20   0.9574468   0.3962264
## 4       0.25   0.9148936   0.4716981
## 5       0.30   0.8936170   0.5471698
## 6       0.35   0.8617021   0.6132075
## 7       0.40   0.8191489   0.6698113
## 8       0.45   0.7553191   0.7169811
## 9       0.50   0.6914894   0.7830189
## 10      0.55   0.6276596   0.8018868
## 11      0.60   0.5744681   0.8584906
## 12      0.65   0.5212766   0.9150943
## 13      0.70   0.4680851   0.9339623
## 14      0.75   0.4148936   0.9339623
## 15      0.80   0.3297872   0.9528302
## 16      0.85   0.2021277   0.9716981
## 17      0.90   0.1276596   0.9716981

Cut-of optimal bisa dipilih berdasarkan: Maksimum dari Sensitivity + Specifcity

8. Catatan ROC cocok saat proporsi kelas seimbang. Untuk data dengan kelas tidak seimbang, precision-recall curve bisa lebih informatif

Precision-Recall Curve (PR Curve)

1. Definisi

Precision (Presisi): Proporsi prediksi positif yang benar-benar positif

\(Precision= \frac{TP}{TP + FP}\)
Recall (Sensitivitas): Proporsi kasus positif yang berhasil diprediksi positif

\(Recall= \frac{TP}{TP + FN}\)

2. Interpretasi

PR Curve menggambarkan trade-off antara Precision dan Recall pada berbagai nilai threshold.
Tujuan ideal: Precision tinggi dan Recall tinggi.
Namun biasanya, peningkatan recall menyebabkan penurunan precision, dan sebaliknya.

3. Area Under PR Curve (AUPRC)

Luas area di bawah kurva PR disebut AUPRC.
AUPRC mendekati 1 → model sangat baik.
Baseline AUPRC = prevalensi kelas positif di data. (Misalnya jika hanya 10% data positif, baseline PR-nya 0.10.)

4. PR Curve vs ROC Curve

Aspek	ROC Curve	Precision-Recall Curve
Fokus	Semua kelas	Kelas positif saja
Cocok untuk	Data seimbang	Data tidak seimbang
Sumbu Y	Sensitivitas / True Positive Rate	Precision
Sumbu X	1 - Spesifisitas / False Positive Rate	Recall

Kapan Gunakan PR Curve?

Gunakan PR Curve jika:

Dataset tidak seimbang (misal: 5% positif, 95% negatif)
Fokus utama adalah deteksi kasus positif (misalnya deteksi penyakit langka, penipuan, dsb.)

5. Visualisasi PR Curve di R

library(PRROC)

## Warning: package 'PRROC' was built under R version 4.3.3

## Loading required package: rlang

pr <- pr.curve(scores.class0 = pred[df$y == 1],
scores.class1 = pred[df$y == 0],
curve = TRUE)
plot(pr)

6. Catatan • PR Curve sangat informatif untuk aplikasi seperti deteksi penipuan atau diagnosis penyakit langka. • Gunakan PR Curve saat: – Kelas positif jauh lebih jarang daripada kelas negatif – Tujuan aplikasi lebih mementingkan presisi terhadap kelas minoritas

Pseudo R-Squared pada Regresi Logistik

Perhitungan Manual R-Squared

logL0 <- logLik(null_model)
logLM <- logLik(model3)
L0 <- exp(logL0)
LM <- exp(logLM)
n <- nobs(model3)
cox_snell <- 1 - (L0 / LM)^(2 / n)
mcfadden <- 1 - (as.numeric(logLM) / as.numeric(logL0))
r2 <- data.frame(
R2_Cox_Snell = cox_snell,
R2_McFadden = mcfadden
)
r2

Perhitungan Otomatis menggunakan Package

if (!require(pscl)) install.packages("pscl"); library(pscl)

## Loading required package: pscl

## Warning: package 'pscl' was built under R version 4.3.3

## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.

pR2(model3)

## fitting null model for pseudo-r2

##          llh      llhNull           G2     McFadden         r2ML         r2CU 
## -102.7365107 -138.2692198   71.0654181    0.2569821    0.2990559    0.3992214

Menggunakan DescTools

if (!require(DescTools)) install.packages("DescTools"); library(DescTools)
PseudoR2(model3, which = "all")

##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.2569821       0.2280530       0.2990559       0.3992214       0.2621707 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.4517796       0.3124594       0.4306149       0.3129335     213.4730215 
##             BIC          logLik         logLik0              G2 
##     226.6662910    -102.7365107    -138.2692198      71.0654181

Interpretasi • Nilai R^2 mendekati 1 berarti model memiliki kekuatan prediktif yang baik. • McFadden R² > 0.2 sering dianggap sebagai model dengan kecocokan yang baik. • Cox & Snell R² lebih konservatif dan tidak pernah mencapai 1 penuh.

Distribusi Multinomial

Definisi

Distribusi multinomial adalah generalisasi dari distribusi binomial ketika hasil dari suatu percobaan dapat jatuh ke dalam lebih dari dua kategori (misalnya, kategori 1 sampai kategori k). Distribusi ini menjelaskan probabilitas dari sejumlah percobaan (n), masing-masing menghasilkan salah satu dari k kategori dengan probabilitas tetap.

\[ P(X_1=x_1,...,X_k = x_k) = \frac{n!}{x_1!x_2!...x_k!}p_1^{x_1}p_2^{x_2}...p_k^{x_k} \]

Studi Kasus

Sebuah hotel ingin memahami preferensi tamu terhadap tiga jenis kamar yang tersedia:

Kamar Standar
Kamar Deluxe
Kamar Suite

Manajemen memiliki data historis bahwa secara umum, proporsi tamu yang memilih masing-masing kamar adalah:

\(p_1=0.5\) untuk Standar
\(p_2=0.3\) untuk Deluxe
\(p_3=0.2\) untuk Suite

Pada suatu hari, ada 10 tamu baru yang memesan kamar, dan hasilnya:

5 tamu memilih Standar
3 tamu memilih Deluxe
2 tamu memilih Suite

Berapa probabilitas bahwa dalam satu hari, terjadi pola pemesanan kamar seperti ini (5 Standar, 3 Deluxe, 2 Suite)?

Perhitungan Manual menggunakan R

n <- 10
x <- c(5, 3, 2)       
p <- c(0.5, 0.3, 0.2) 

faktorial_total <- factorial(n)
faktorial_x <- prod(factorial(x))
koefisien <- faktorial_total / faktorial_x
peluang <- koefisien * prod(p^x)
peluang

## [1] 0.08505

Probabilitas bahwa 5 orang memilih Kamar Standar, 3 Kamar Deluxe, dan 2 Kamar Suite (dengan proporsi preferensi 0.5, 0.3, dan 0.2) adalah 0.08505. Distribusi multinomial digunakan untuk menghitung peluang dalam percobaan dengan beberapa kategori hasil. Rumus dasarnya merupakan generalisasi dari binomial untuk lebih dari dua kategori.

Multinomial Regresi Logistik

Regresi logistik multinomial digunakan untuk memodelkan hubungan antara satu variabel respon kategorik dengan lebih dari dua kategori dan satu atau lebih variabel prediktor.

Misalkan Y memiliki K kategori:

\(Y∈{1,2,…,K}\)

Kita pilih salah satu kategori sebagai referensi (baseline), misalnya kategori K. Maka, untuk setiap kategori \(j=1,2,…,K−1\) , kita bangun model:

\(log⁡(\frac{P(Y=j)}{P(Y=K)})=βj0+βj1X1+βj2X2+⋯+βjpXp\)

\(P(Y=j)\) adalah probabilitas bahwa respons berada pada kategori ke-j
\(P(Y=k)\) adalah probabilitas kategori referensi

Baseline-category logit model

Baseline-category logit model digunakan dalam regresi logistik multinomial (nominal), untuk memodelkan probabilitas kategori respon dengan membandingkan setiap kategori terhadap satu kategori referensi (baseline).

Rumus Umum

Untuk K kategori, dan memilih satu kategori (misal ke-K) sebagai baseline:

\(log⁡(\frac{πj}{πK})=βj0+βj1X,j=1,…,K−1\)

\(π_j\): probabilitas respon berada di kategori j
\(π_K\) : probabilitas respon berada di kategori baseline K
X: prediktor
Jumlah fungsi logit: K-1

Estimasi Parameter

Estimasi dilakukan menggunakan metode Maximum Likelihood (ML), yang bertujuan untuk menemukan parameter \(\theta\) yang memaksimalkan log-likelihood function berdasarkan data pengamatan.

Log-likelihood function:

\(ℓ(θ)=∑y_{ij} \log⁡(\pi_{ij})\)

Studi Kasus

Sebuah dinas perhubungan kota ingin memahami faktor-faktor yang memengaruhi preferensi warga dalam memilih moda transportasi utama yang digunakan untuk pergi bekerja atau beraktivitas.

Dinas ini melakukan survei terhadap 200 responden dan mencatat informasi berikut:

Transportasi: Moda transportasi yang dipilih (Mobil, Motor, Transportasi Umum)
Pendapatan: Pendapatan bulanan dalam juta rupiah
Usia: Usia responden
Jarak: Jarak tempuh harian dari rumah ke tempat kerja/sekolah (dalam km)

Tujuan Analisis: Mengetahui bagaimana pendapatan, usia, dan jarak tempuh memengaruhi pilihan transportasi utama warga kota.

Simulasi Data

set.seed(456)

# Jumlah responden
n <- 200

# Variabel independen
Pendapatan <- round(rnorm(n, mean = 6, sd = 2), 1) # Pendapatan 6 juta rerata
Usia <- round(rnorm(n, mean = 35, sd = 10))
Jarak <- round(abs(rnorm(n, mean = 10, sd = 5)), 1) # Jarak dalam km

# Simulasi preferensi berdasarkan pendapatan dan jarak
Transportasi <- sapply(1:n, function(i) {
  if (Pendapatan[i] >= 8) {
    sample(c("Mobil", "Motor", "Transportasi Umum"), 1, prob = c(0.6, 0.3, 0.1))
  } else if (Pendapatan[i] >= 4) {
    sample(c("Mobil", "Motor", "Transportasi Umum"), 1, prob = c(0.3, 0.5, 0.2))
  } else {
    sample(c("Mobil", "Motor", "Transportasi Umum"), 1, prob = c(0.1, 0.4, 0.5))
  }
})

# Buat data.frame
df <- data.frame(
  Transportasi = factor(Transportasi),
  Pendapatan = Pendapatan,
  Usia = Usia,
  Jarak = Jarak
)

# Set baseline kategori
df$Transportasi <- relevel(df$Transportasi, ref = "Transportasi Umum")

# Tampilkan sebagian data
head(df)

Estimasi Model

library(nnet) 
model_mnlogit <- multinom(Transportasi ~ Pendapatan + Usia + Jarak, data = df)

## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 194.979797
## final  value 194.957424 
## converged

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) Pendapatan   Usia  Jarak
## Mobil      0.0001     0.0000 0.0631 0.7250
## Motor      0.0642     0.0044 0.1770 0.7749

Interpretasi :

Niali p-value < 0.05 berarti bahwa variabel tersebut signifikan terhadap preferensi Transportasi. Pada kasus ini Pendapatan memengaruhi preferensi transportasi.

Prediksi dan Validasi

df$Predicted <- predict(model_mnlogit)
table(Predicted = df$Predicted, Actual = df$Transportasi)

##                    Actual
## Predicted           Transportasi Umum Mobil Motor
##   Transportasi Umum                14     2     6
##   Mobil                             7    34    16
##   Motor                            26    30    65

Visualisasi Plot

ggplot(df, aes(x = Pendapatan, y = Jarak, fill = Predicted)) +
  geom_tile(alpha = 0.6) +
  geom_point(data = df, aes(x = Pendapatan, y = Jarak, color = Transportasi), size = 2, alpha = 0.7) +
  scale_fill_manual(values = c("Mobil" = "#E41A1C", "Motor" = "#377EB8", "Transportasi Umum" = "#4DAF4A")) +
  scale_color_manual(values = c("Mobil" = "#E41A1C", "Motor" = "#377EB8", "Transportasi Umum" = "#4DAF4A")) +
  labs(title = "Prediksi Moda Transportasi Berdasarkan Pendapatan dan Jarak",
       x = "Pendapatan (juta Rupiah)", y = "Jarak (km)",
       fill = "Prediksi", color = "Data Aktual") +
  theme_minimal()

Kesimpulan

Model regresi logistik multinomial berhasil digunakan untuk:

Menganalisis hubungan antara atribut karyawan dan preferensi transportasi
Mengetahui faktor signifkan yang memengaruhi pilihan transportasi
Memungkinkan prediksi transportasi oleh karyawan baru berdasarkan karakteristiknya

Regresi logistik ordinal

Regresi logistik ordinal digunakan ketika variabel respon Y bersifat ordinal (memiliki urutan), misalnya tingkat kepuasan: Rendah, Sedang, Tinggi.

Model ini berbeda dengan:

Regresi logistik biner: hanya 2 kategori
Regresi logistik multinomial: kategori > 2 tetapi tidak berurutan

Konsep Cumulative Logit Model

Model ini digunakan untuk data ordinal—yaitu data kategorik yang memiliki urutan tetapi tidak memiliki jarak yang pasti antar kategori (contoh: tingkat kepuasan, tingkat pendidikan, dsb).

Bentuk Model:

\(log⁡(\frac{P(Y≤j)}{P(Y>j)})=αj+βX\)

Y: variabel respon ordinal dengan K kategori
j=1,2,…,K−1: titik batas antar kategori
αj: intercept khusus untuk cutoff ke-j
β: koefisien regresi (sama untuk semua model → proportional odds.
X: variabel prediktor

Model ini membentuk K−1 model logit kumulatif

Interpretasi Koefisien

Jika β>0:

-   Semakin besar nilai X, semakin besar peluang untuk berada di **kategori lebih rendah** (misalnya dari "Sangat Baik" ke "Baik"

Jika β<0:
- Semakin besar nilai X, semakin besar peluang untuk berada di kategori lebih tinggi (dalam urutan).

Odds Ratio (OR)

Untuk interpretasi mudah:

\[ OR=e^\beta \]

Jika OR > 1, kenaikan 1 unit di X meningkatkan peluang berada di kategori rendah atau sama.
Jika OR < 1, kenaikan 1 unit di X meningkatkan peluang ke kategori tinggi.

Studi Kasus

X1 : Usia

Y = Pendapat mengenai Kecepatan Pelayanan Rumah Sakit (1= Sangat Lama, 2 = Lama, 3 = Sedang, 4 = Cepat, 5 = Sangat Cepat)

Tujuannya adalah: Mengetahui bagaimana Usia memengaruhi Pendapat mengenai Kecepatan Pelayanan Rumah Sakit

Umur <- c(35,54,32,33,26,23,30,39,27,54,42,21,57,57,42,53,18,38,51,32,28,19,36,48,29,50,50,39,45,32,57,51,45,21,53,39,31,52,46,29,49,49,48,49,20,37,64,24,23,25,32,28,47,35,34,29,47,39,40)  

ResponAsli <- c(4,4,3,4,4,3,4,5,4,3,5,5,3,5,2,4,4,3,1,4,3,4,3,5,5,3,5,1,4,3,4,4,4,3,4,3,4,4,4,3,3,1,3,3,4,2,3,3,3,4,4,3,5,5,4,5,4,3,3) 

# Buat data.frame 
df <- data.frame(Umur = Umur,ResponAsli = ResponAsli)  
# Buat variabel kategori baru: ResponCepat 
df$ResponCepat <- with(df, case_when(
  ResponAsli == 1 ~ "Sangat Lambat",   
  ResponAsli == 2 ~ "Lambat",   
  ResponAsli == 3 ~ "Sedang",   
  ResponAsli == 4 ~ "Cepat",  
  ResponAsli == 5 ~ "Sangat Cepat",   TRUE ~ NA_character_ ))  
# Konversi ResponCepat ke faktor dan set 
baseline = "Sangat Lambat" 
df$ResponCepat <- factor(df$ResponCepat, levels = c("Sangat Lambat", "Lambat", "Sedang", "Cepat", "Sangat Cepat")) 
# Tampilkan 6 baris pertama 
head(df)

Estimasi Model Ordinal

model_ord <- polr(ResponCepat ~ Umur, data = df, Hess = TRUE)
summary(model_ord)

## Call:
## polr(formula = ResponCepat ~ Umur, data = df, Hess = TRUE)
## 
## Coefficients:
##          Value Std. Error t value
## Umur -0.006789    0.01997 -0.3399
## 
## Intercepts:
##                      Value   Std. Error t value
## Sangat Lambat|Lambat -3.1898  0.9771    -3.2646
## Lambat|Sedang        -2.6413  0.9034    -2.9236
## Sedang|Cepat         -0.4969  0.8047    -0.6175
## Cepat|Sangat Cepat    1.3312  0.8322     1.5997
## 
## Residual Deviance: 153.5152 
## AIC: 163.5152

Nilai P-Value dan Interpretasi

(ctable <- coef(summary(model_ord)))

##                             Value Std. Error    t value
## Umur                 -0.006788812 0.01997067 -0.3399392
## Sangat Lambat|Lambat -3.189755048 0.97706672 -3.2646236
## Lambat|Sedang        -2.641296099 0.90342856 -2.9236358
## Sedang|Cepat         -0.496895403 0.80473416 -0.6174653
## Cepat|Sangat Cepat    1.331211251 0.83216476  1.5996967

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
## Umur                 -0.006788812 0.01997067 -0.3399392  0.7339
## Sangat Lambat|Lambat -3.189755048 0.97706672 -3.2646236  0.0011
## Lambat|Sedang        -2.641296099 0.90342856 -2.9236358  0.0035
## Sedang|Cepat         -0.496895403 0.80473416 -0.6174653  0.5369
## Cepat|Sangat Cepat    1.331211251 0.83216476  1.5996967  0.1097

Interpretasi :

Niali p-value < 0.05 berarti bahwa variabel tersebut signifikan terhadap respon kecepatan pelayanan di RSHS. Pada kasus ini Usia tidak memengaruhi kecepatan pelayanan. Sedangkan threshold Sangat Lambat | Lambat dan Lambat | Sedang berpengaruh signidikan

Prediksi Probablitias

newdata <- data.frame(Umur = 5:9)
predict(model_ord, newdata = newdata, type = "probs")

##   Sangat Lambat     Lambat    Sedang     Cepat Sangat Cepat
## 1    0.04086292 0.02780383 0.3176192 0.4103104    0.2034037
## 2    0.04112982 0.02797236 0.3187944 0.4097975    0.2023059
## 3    0.04139840 0.02814177 0.3199695 0.4092777    0.2012126
## 4    0.04166865 0.02831207 0.3211445 0.4087511    0.2001236
## 5    0.04194059 0.02848326 0.3223193 0.4082177    0.1990391

Goodness of Fit dan Proportional Odds

Model cumulative logit mengasumsikan efek prediktor sama untuk setiap cutof. Jika tidak, pertimbangkan model non-proportional odds seperti generalized ordinal model.

Alternatif Model Ordinal

Selain cumulative logit, model ordinal lainnya: • Adjacent-category logit • Continuation-ratio (sequential) logit Model tersebut dapat digunakan saat asumsi proportional odds tidak terpenuhi.

Kesimpulan

Regresi ordinal efektif untuk respon berurutan.
Model cumulative logit menginterpretasikan efek dalam bentuk log-odds kumulatif.
Implementasi di R dengan fungsi polr() dari package MASS.

Jika diperlukan validasi lanjut, bisa digunakan uji devians atau likelihood ratio test.

Asumsi Paralelisme dalam Regresi Logistik Ordinal

Asumsi paralelisme atau proportional odds dalam model cumulative logit menyatakan bahwa:

Efek prediktor (koefisien β) konstan untuk semua batas kumulatif kategori.

Dengan kata lain, hanya intersep (\(\alpha_j\)) yang berubah untuk setiap batas kategori, sedangkan kemiringan (β) tetap.

Bentuk Model Matematis

Untuk variabel respon ordinal Y dengan K kategori, model cumulative logit adalah:

\(log⁡(frac{P(Y≤j)}{P(Y>j)})=α_j+β_X\)

Untuk j=1,2,…,K−1

Interpretasi

Jika β>0: prediktor meningkatkan peluang berada di kategori lebih rendah.
Jika β<0 : prediktor meningkatkan peluang berada di kategori lebih tinggi.
Odds ratio: exp(β)

Visualisasi

Dalam asumsi ini, grafik logit kumulatif terhadap prediktor akan tampak sejajar (paralel), menunjukkan efek prediktor konsisten di seluruh tingkat kategori.

Konsekuensi Jika Asumsi Dilanggar

Jika asumsi tidak terpenuhi:

Efek prediktor berbeda untuk tiap batas kategori.
Interpretasi menjadi tidak valid.

Solusi:

Gunakan model alternatif seperti:

-    **Generalized Ordinal Logistic Regression**

-   **Partial Proportional Odds Model**

Pengujian Asumsi Paralelisme

1. Likelihood Ratio Test (Manual)

Bandingkan:

Model dengan asumsi proportional odds (polr)
Model non-proportional (misal multinom() dari nnet)

2. Brant Test (Paket brant)

Jika p-value < 0.05 → asumsi paralelisme ditolak

Kesimpulan

Asumsi paralelisme sangat penting untuk validitas model cumulative logit.
Jika terpenuhi, model lebih sederhana dan mudah diinterpretasikan
Jika tidak, gunakan model ordinal alternatif agar interpretasi tetap sah.

Log Linear Model

Dalam berbagai bidang ilmu terapan seperti kesehatan masyarakat, ilmu sosial, pemasaran, dan epidemiologi, kita sering kali berhadapan dengan data kategorik, yaitu data yang diklasifikasikan ke dalam kelompok atau kategori tertentu. Contohnya meliputi:

Jenis kelamin (laki-laki/perempuan),
Status merokok (pernah/tidak pernah),
Tingkat pendidikan (SMA/S1/S2),
Status kesehatan (sakit/sehat),
Preferensi produk (merek A/B/C).

Data semacam ini tidak bisa dianalisis dengan metode statistik parametrik biasa seperti regresi linear karena asumsi yang mendasarinya—seperti normalitas residual dan hubungan linier antar variabel—tidak terpenuhi pada data kategorik. Oleh karena itu, dibutuhkan pendekatan khusus untuk menganalisis hubungan antar variabel kategorik. Untuk menganalisis data seperti ini, digunakan:

Tabel Kontingensi

Tabel kontingensi menyajikan jumlah frekuensi dari kombinasi kategori antar variabel.

table_data <- matrix(c(30, 20, 50, 70), nrow=2,
dimnames = list(Pendidikan = c("SMP", "SMA"),
KecepatanPelayanan = c("Ya", "Tidak")))
table_data

##           KecepatanPelayanan
## Pendidikan Ya Tidak
##        SMP 30    50
##        SMA 20    70

Tabel kontingensi bersifat deskriptif dan tidak melibatkan pemodelan probabilitas.

Model Log-linear

Model loglinear memodelkan logaritma dari ekspektasi frekuensi sel dalam tabel kontingensi. \[log(\mu_{ij})=\mu+\lambda_i^A+\lambda_j^B+\lambda_{ij}^{AB }\] • Cocok untuk menyelidiki asosiasi dan independensi antar variabel. • Tidak membedakan antara variabel respon dan penjelas. • Umumnya digunakan pada tabel multi-dimensi (2 arah, 3 arah, dst).
```
library(MASS)
loglm(~ Pendidikan * KecepatanPelayanan, data = table_data)
```
```
## Call:
## loglm(formula = ~Pendidikan * KecepatanPelayanan, data = table_data)
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1
```
Model Regresi Logistik

Model regresi logistik biner:
$\log(\frac{p}{1-p})=\beta_0+\beta_1x$
• Digunakan jika ada variabel dependen kategorik (biasanya biner).
• Bertujuan untuk memprediksi probabilitas suatu outcome.
• Umumnya digunakan dalam studi observasional atau eksperimental.


``` r
data_glm <- data.frame(
KecepatanPelayanan = c(1, 0, 1, 0),
Pendidikan = factor(c("SMP", "SMP", "SMA", "SMA")),
Frek = c(30, 20, 50, 70)
)
model_logit <- glm(KecepatanPelayanan ~ Pendidikan, weights = Frek, family = binomial, data = data_glm)
summary(model_logit)
```

```
## 
## Call:
## glm(formula = KecepatanPelayanan ~ Pendidikan, family = binomial, 
##     data = data_glm, weights = Frek)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    -0.3365     0.1852  -1.817   0.0692 .
## PendidikanSMP   0.7419     0.3430   2.163   0.0305 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 235.08  on 3  degrees of freedom
## Residual deviance: 230.31  on 2  degrees of freedom
## AIC: 234.31
## 
## Number of Fisher Scoring iterations: 4
```

Pendekatan	Tujuan	Variabel Dependen	Fokus Analisis	Distribusi
Tabel Kontingensi	Deskriptif eksploratif	Tidak	Hubungan antar kategori	Tidak digunakan
Model Log-linear	Eksplorasi asosiasi	Tidak	Interaksi/simetris antar variabel	Poisson
Regresi Logistik	Prediktif	Ya	Prediksi probabilitas outcome	Binomial

Tabel Kontingensi dan Model Loglinear

Tabel kontingensi menyajikan jumlah frekuensi dari kombinasi kategori antar variabel.

table_data <- matrix(c(30, 20, 50, 70), nrow=2,
dimnames = list(Pendidikan = c("SMP", "SMA"),
KecepatanPelayanan = c("Ya", "Tidak")))
table_data

##           KecepatanPelayanan
## Pendidikan Ya Tidak
##        SMP 30    50
##        SMA 20    70

Model loglinear memodelkan logaritma dari ekspektasi frekuensi sel dalam tabel kontingensi. \[log(\mu_{ij})=\mu+\lambda_i^A+\lambda_j^B+\lambda_{ij}^{AB }\]

Model Saturated

Model saturated atau model penuh menyertakan seluruh efek utama dan interaksi: • Cocok sempurna terhadap data • Tidak mengasumsikan independensi antar variabel Contoh formulasi untuk tabel 2x2:

model_saturated <- loglm(~ KecepatanPelayanan * Pendidikan, data = table_data)
summary(model_saturated)

## Formula:
## ~KecepatanPelayanan * Pendidikan
## attr(,"variables")
## list(KecepatanPelayanan, Pendidikan)
## attr(,"factors")
##                    KecepatanPelayanan Pendidikan KecepatanPelayanan:Pendidikan
## KecepatanPelayanan                  1          0                             1
## Pendidikan                          0          1                             1
## attr(,"term.labels")
## [1] "KecepatanPelayanan"            "Pendidikan"                   
## [3] "KecepatanPelayanan:Pendidikan"
## attr(,"order")
## [1] 1 1 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Model Independent

Model independen mengasumsikan bahwa tidak ada interaksi antara variabel: \[\log(\mu_{ij})=\mu+\lambda_i^X+\lambda_j^Y\] Model ini menguji hipotesis bahwa variabel X dan Y saling independen.

model_indep <- loglm(~ KecepatanPelayanan + Pendidikan, data = table_data)
summary(model_indep)

## Formula:
## ~KecepatanPelayanan + Pendidikan
## attr(,"variables")
## list(KecepatanPelayanan, Pendidikan)
## attr(,"factors")
##                    KecepatanPelayanan Pendidikan
## KecepatanPelayanan                  1          0
## Pendidikan                          0          1
## attr(,"term.labels")
## [1] "KecepatanPelayanan" "Pendidikan"        
## 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 4.773916  1 0.02889403
## Pearson          4.761574  1 0.02910192

Odds Ratio dan Interpretasi

Odds ratio untuk tabel 2x2: \[OR=\frac{n_{11}n_{22}}{n_{12}n_{21}}\] Interpretasi nilai OR: • OR = 1: Tidak ada asosiasi • OR > 1: Asosiasi positif • OR < 1: Asosiasi negatif

Estimasi Parameter

Dalam model saturated: • Estimasi dilakukan dengan pembatasan seperti sum-to-zero • Estimasi parameter dilakukan dengan iterative proportional ftting (IPF)

# Estimasi odds ratio dan log-odds
logOR <- log((table_data[1,1] * table_data[2,2]) / (table_data[1,2] * table_data[2,1]))
logOR

## [1] 0.7419373

Model Lebih Sederhana dan Perbandingan Model

Perbandingan antar model dilakukan dengan menggunakan statistik deviance (G²) atau likelihood ratio test.

anova(model_indep, model_saturated)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~KecepatanPelayanan + Pendidikan 
## Model 2:
##  ~KecepatanPelayanan * Pendidikan 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   4.773916  1                                    
## Model 2   0.000000  0   4.773916         1        0.02889
## Saturated 0.000000  0   0.000000         0        1.00000

Data Survey

Suatu Penelitian ingin meneliti seberapa besar pengaruh Tingkat Pendidikan terhadap persepsi lamanya pelayanan Rumah Sakit Hasan Sadikin. Survei ini dilakukan pada tanggal 3 - 4 Juni 2025 dan dilakukan pada Ruang Tunggu IGD, ICU, Loket Pengambilan Hasil Lab, dan Radiologi. Terdapat 33 Responden yang diwawancarai.

data_survey <- matrix(c(10,0,
13,2,
3,3,
2,0),
nrow = 4, byrow = TRUE,
dimnames = list(Pendidikan = c("Tinggi", "SMA", "SMP","SD"),
KecepatanPelayanan = c("Cepat", "Tidak Cepat")))
ftable(data_survey)

##            KecepatanPelayanan Cepat Tidak Cepat
## Pendidikan                                     
## Tinggi                           10           0
## SMA                              13           2
## SMP                               3           3
## SD                                2           0

loglm(~ Pendidikan + KecepatanPelayanan, data = data_survey)

## Call:
## loglm(formula = ~Pendidikan + KecepatanPelayanan, data = data_survey)
## 
## Statistics:
##                       X^2 df   P(> X^2)
## Likelihood Ratio 7.973667  3 0.04655907
## Pearson          7.849286  3 0.04923140

Interpretasi : Ada hubungan yang signifikan antara tingkat pendidikan dan persepsi terhadap kecepatan pelayanan. Karena p-va;ue < \[\alpha=0.05\]

Model Log-Linear Tiga Arah

Dalam dunia nyata, banyak fenomena melibatkan lebih dari dua variabel kategorik. Misalnya, dalam studi sosial, pendidikan, dan kesehatan, sering kali kita tertarik pada hubungan kompleks antara tiga atau lebih faktor, seperti jenis kelamin, tingkat pendidikan, dan status pekerjaan, atau dalam konteks medis seperti jenis terapi, usia pasien, dan hasil pengobatan. Untuk menganalisis hubungan ini, model log-linear menjadi alat yang sangat efektif.

Model log-linear tiga arah adalah perluasan dari model dua arah yang digunakan untuk mengkaji struktur hubungan di antara tiga variabel kategorik secara simultan. Model ini tidak menganggap adanya variabel dependen, sehingga semua variabel dianggap setara (simetris), berbeda dari model regresi logistik yang fokus pada prediksi satu variabel dependen.

Tujuan utama model log-linear adalah untuk mengidentifikasi apakah terdapat asosiasi antara ketiga variabel, dan sejauh mana interaksi antar variabel tersebut memengaruhi distribusi frekuensi dalam tabel kontingensi.

Model Log-Linear untuk tabel Tiga Arah

1. Pengertian

Model log-linear tiga arah digunakan untuk menganalisis hubungan antara tiga variabel kategorik dalam tabel kontingensi tiga dimensi. Berbeda dengan regresi logistik, model ini tidak membedakan variabel dependen dan independen, karena semua variabel dianggap simetris.

2. Tujuan

Mengidentifikasi asosiasi dan interaksi antar tiga variabel kategorik.
Menentukan apakah interaksi dua arah atau tiga arah diperlukan untuk menjelaskan struktur data.

3. Struktur Model

Model dinyatakan sebagai:

\(log⁡(μ_{ijk})=λ+λ_i^X+λ_j^Y+λ_k^Z+λ_{ij}^{XY}+λ_{ik}^{XZ}+λ_{jk}^{YZ}+λ_{ijk}^{XYZ}\)

4. Jenis Model

Model	Ciri Utama
Saturated	Semua efek utama, dua arah, dan tiga arah
Homogeneous Association	Hanya efek utama + semua interaksi dua arah
Conditional Independence	Interaksi dua arah tertentu saja
Joint Independence	Hanya satu interaksi dua arah
Mutual Independence	Hanya efek utama (tanpa interaksi)

5. Manfaat

Menjelaskan struktur hubungan antar tiga variabel.
Menguji apakah interaksi dua atau tiga arah signifikan.
Cocok untuk studi sosial, pendidikan, kesehatan, dan pasar.

6. Pengujian Model

Gunakan Likelihood Ratio Test atau Chi-Square untuk menilai kecocokan model.
Bandingkan model sederhana dan kompleks untuk melihat pengaruh interaksi

Menilai hubungan antara:

Jenis Kelamin (X)
Pendidikan (Y)
Status Pekerjaan (Z)

Apakah hubungan Pendidikan - Pekerjaan berbeda antar jenis kelamin?

Pengujian Interaksi dalam Model Log-Linear Tiga Arah

Dalam analisis model log-linear tiga arah, pengujian interaksi dilakukan untuk mengetahui ada atau tidaknya interaksi antar variabel. Pengujian ini dilakukan secara bertahap, dimulai dari tingkat interaksi tertinggi ke yang lebih rendah. Untuk model log-linear dengan tiga peubah (X, Y, dan Z), tahapan pengujian meliputi:

Pengujian Interaksi Tiga Arah (XYZ):
- Bandingkan model saturated dengan model homogenous.
Pengujian Interaksi Dua Arah (XY, XZ, YZ):
- Bandingkan model homogenous dengan model conditional.
- Bandingkan model conditional dengan model joint independence.
- Bandingkan model joint independence dengan model tanpa interaksi.

Setiap tahapan pengujian dilakukan untuk menilai kecocokan model dan menentukan struktur interaksi mana yang paling sesuai dengan data yang diamati.

Studi Kasus

Suatu Penelitian ingin meneliti seberapa besar pengaruh Jenis Kelamin dan Tingkat Pendidikan terhadap persepsi lamanya pelayanan Rumah Sakit Hasan Sadikin. Survei ini dilakukan pada tanggal 3 - 4 Juni 2025 dan dilakukan pada Ruang Tunggu IGD, ICU, Loket Pengambilan Hasil Lab, dan Radiologi.

X : Jenis Kelamin (Laki/Perempuan)

Z : Tingkat Pendidikan (SD, SMP, SMA, Tinggi)

Y = Kecepatan Pelayanan Rumah Sakit (0 = Lama / Tidak Cepat, 1 = Cepat )

Data Survei
Pendidikan	Jenis Kelamin	Cepat	Tidak Cepat
Tinggi	Laki-Laki	6	0
	Perempuan	4	0
SMA	Laki-Laki	5	2
	Perempuan	8	0
SMP	Laki-Laki	2	1
	Perempuan	3	0
SD	Laki-Laki	1	0
	Perempuan	1	0

Analisis Log-Linear untuk Tabel Tiga Arah

Package yang Digunakan

library("epitools")
library("DescTools")
library("lawstat")

## Warning: package 'lawstat' was built under R version 4.3.3

# Input data sesuai tabel praktikum
z.edu <- factor(rep(c("Tinggi", "SMA", "SMP","SD"), each = 4))
x.sex  <- factor(rep(c("Laki", "Perempuan"), each = 2, times = 4))
y.fav  <- factor(rep(c("Cepat", "Tidak Cepat"), times = 8))
counts <- c(6,0,4,0,5,2,8,0,2,0,3,1,1,0,1,0)

data <- data.frame(
  Pendidikan = z.edu,
  Jenis_Kelamin   = x.sex,
  Sikap           = y.fav,
  Frekuensi       = counts

)
data

Membentuk Tabel Kontingensi 3 Arah

table3d <- xtabs(Frekuensi ~ Pendidikan + Jenis_Kelamin + Sikap, data = data)
ftable(table3d)

##                          Sikap Cepat Tidak Cepat
## Pendidikan Jenis_Kelamin                        
## SD         Laki                    1           0
##            Perempuan               1           0
## SMA        Laki                    5           2
##            Perempuan               8           0
## SMP        Laki                    2           0
##            Perempuan               3           1
## Tinggi     Laki                    6           0
##            Perempuan               4           0

UJI MODEL INTERAKSI TIGA ARAH (SATURATED VS HOMOGENOUS)

Penentuan Kategori Referensi

##=============================##
# Penentuan kategori reference
##=============================##
x.sex <- relevel(x.sex, ref = "Perempuan")
y.fav <- relevel(y.fav, ref = "Tidak Cepat")
z.fund <- relevel(z.edu, ref = "SD")

Model Saturated

Model log-linear 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(counts ~ x.sex + y.fav + z.edu +
x.sex*y.fav + x.sex*z.edu + y.fav*z.edu +
x.sex*y.fav*z.edu,
family = poisson(link = "log"))
summary(model_saturated)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.edu + x.sex * y.fav + 
##     x.sex * z.edu + y.fav * z.edu + x.sex * y.fav * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                                    Estimate Std. Error z value Pr(>|z|)
## (Intercept)                      -2.430e+01  1.148e+05       0        1
## x.sexLaki                         6.187e-15  1.624e+05       0        1
## y.favCepat                        2.430e+01  1.148e+05       0        1
## z.eduSMA                          1.401e-09  1.624e+05       0        1
## z.eduSMP                          2.430e+01  1.148e+05       0        1
## z.eduTinggi                       7.373e-10  1.624e+05       0        1
## x.sexLaki:y.favCepat              3.693e-15  1.624e+05       0        1
## x.sexLaki:z.eduSMA                2.500e+01  1.989e+05       0        1
## x.sexLaki:z.eduSMP               -2.430e+01  1.989e+05       0        1
## x.sexLaki:z.eduTinggi            -6.196e-10  2.297e+05       0        1
## y.favCepat:z.eduSMA               2.079e+00  1.624e+05       0        1
## y.favCepat:z.eduSMP              -2.320e+01  1.148e+05       0        1
## y.favCepat:z.eduTinggi            1.386e+00  1.624e+05       0        1
## x.sexLaki:y.favCepat:z.eduSMA    -2.547e+01  1.989e+05       0        1
## x.sexLaki:y.favCepat:z.eduSMP     2.390e+01  1.989e+05       0        1
## x.sexLaki:y.favCepat:z.eduTinggi  4.055e-01  2.297e+05       0        1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4.6315e+01  on 15  degrees of freedom
## Residual deviance: 3.3474e-10  on  0  degrees of freedom
## AIC: 60.561
## 
## Number of Fisher Scoring iterations: 22

exp(model_saturated$coefficients)

##                      (Intercept)                        x.sexLaki 
##                     2.789468e-11                     1.000000e+00 
##                       y.favCepat                         z.eduSMA 
##                     3.584913e+10                     1.000000e+00 
##                         z.eduSMP                      z.eduTinggi 
##                     3.584913e+10                     1.000000e+00 
##             x.sexLaki:y.favCepat               x.sexLaki:z.eduSMA 
##                     1.000000e+00                     7.169826e+10 
##               x.sexLaki:z.eduSMP            x.sexLaki:z.eduTinggi 
##                     2.789468e-11                     1.000000e+00 
##              y.favCepat:z.eduSMA              y.favCepat:z.eduSMP 
##                     8.000000e+00                     8.368404e-11 
##           y.favCepat:z.eduTinggi    x.sexLaki:y.favCepat:z.eduSMA 
##                     4.000000e+00                     8.717088e-12 
##    x.sexLaki:y.favCepat:z.eduSMP x.sexLaki:y.favCepat:z.eduTinggi 
##                     2.389942e+10                     1.500000e+00

Ringkasan Model

Model yang digunakan adalah model log-linear saturated dengan semua efek utama, interaksi dua arah, dan interaksi tiga arah. Model ini memodelkan hubungan antara jenis kelamin (x.sex), sikap terhadap kecepatan pelayanani (y.fav), dan tingkat pendidikan(z.edu) terhadap frekuensi responden.

Hasil Estimasi Koefisien

Tabel Hasil Estimasi Koefisien (Model Saturated)

No	Parameter	Estimate	Std. Error	Pr(>
1	(Intercept)	-24.30	114800.00	1.000
2	x.sexLaki	6.19e-15	162400.00	1.000
3	y.favCepat	24.30	114800.00	1.000
4	z.eduSMA	1.40e-09	162400.00	1.000
5	z.eduSMP	24.30	114800.00	1.000
6	z.eduTinggi	7.37e-10	162400.00	1.000
7	x.sexLaki:y.favCepat	3.69e-15	162400.00	1.000
8	x.sexLaki:z.eduSMA	25.00	198900.00	1.000
9	x.sexLaki:z.eduSMP	-24.30	198900.00	1.000
10	x.sexLaki:z.eduTinggi	-6.20e-10	229700.00	1.000
11	y.favCepat:z.eduSMA	2.08	162400.00	1.000
12	y.favCepat:z.eduSMP	-23.20	114800.00	1.000
13	y.favCepat:z.eduTinggi	1.39	162400.00	1.000
14	x.sexLaki:y.favCepat:z.eduSMA	-25.47	198900.00	1.000
15	x.sexLaki:y.favCepat:z.eduSMP	23.90	198900.00	1.000
16	x.sexLaki:y.favCepat:z.eduTinggi	0.41	229700.00	1.000

Interpretasi koefisien

Semua parameter memilki p-value yang sama yaitu 1 dimana p-value > alpha sehingga tidak ada parameter yang signifikan.

Goodness of Fit

Residual deviance≈0 menandakan model saturated benar-benar fit terhadap data (seluruh variasi data dijelaskan oleh model).
AIC = 60.561 dapat digunakan untuk perbandingan dengan model yang lebih sederhana

Kesimpulan

Model saturated ini sangat fit dengan data, namun semua parameter/interaksi tidak signifikan.
Tidak ditemukan bukti kuat interaksi dua atau tiga arah yang signifikan.
Model yang lebih sederhana (tanpa interaksi tiga arah) perlu dipertimbangkan untuk model final yang lebih parsimonious.

Model Homogenous

Model log-linear homogenous memasukkan semua efek utama dan semua interaksi dua arah, tanpa interaksi tiga arah. Secara matematis, model ini dapat dituliskan sebagai berikut:

\[ \log(\mu_{ijk})=\lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{XZ}+\lambda_{}jk^{YZ} \]

# Homogenous Model
model_homogenous <- glm(counts ~ x.sex + y.fav + z.edu +
              x.sex*y.fav + x.sex*z.edu + y.fav*z.edu,
              family = poisson(link = "log"))
summary(model_homogenous)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.edu + x.sex * y.fav + 
##     x.sex * z.edu + y.fav * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)              -21.0188 10681.1527  -0.002    0.998
## x.sexLaki                  1.2040     1.9463   0.619    0.536
## y.favCepat                21.0188 10681.1526   0.002    0.998
## z.eduSMA                  20.4592 10681.1527   0.002    0.998
## z.eduSMP                  20.1715 10681.1527   0.002    0.998
## z.eduTinggi               -0.2585 14859.0930   0.000    1.000
## x.sexLaki:y.favCepat      -1.2040     1.3372  -0.900    0.368
## x.sexLaki:z.eduSMA        -0.2877     1.5171  -0.190    0.850
## x.sexLaki:z.eduSMP        -0.9163     1.6904  -0.542    0.588
## x.sexLaki:z.eduTinggi      0.4055     1.5546   0.261    0.794
## y.favCepat:z.eduSMA      -18.4539 10681.1526  -0.002    0.999
## y.favCepat:z.eduSMP      -18.8985 10681.1527  -0.002    0.999
## y.favCepat:z.eduTinggi     1.6448 14859.0930   0.000    1.000
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 46.3151  on 15  degrees of freedom
## Residual deviance:  3.4438  on  3  degrees of freedom
## AIC: 58.005
## 
## Number of Fisher Scoring iterations: 18

Uji Hipotesis: Apakah Ada Interaksi Tiga Arah? (Saturated vs Homogenous)

Pengujian ini menggunakan residual deviance dari kedua model (saturated dan homogenous).

Langkah-Langkah Pengujian

1. Hipotesis

H0: Tidak ada interaksi tiga arah (model homogenous sudah cukup)
H1: Ada interaksi tiga arah (model saturated diperlukan)

2. Hitung Selisih Deviance

# Deviance antar model
Deviance.model <- model_homogenous$deviance - model_saturated$deviance
Deviance.model

## [1] 3.443845

3. Hitung Derajat Bebas

# Derajat bebas = db model homogenous - db model saturated
derajat.bebas <- (model_homogenous$df.residual - model_saturated$df.residual)
derajat.bebas

## [1] 3

4. Chi-Square Tabel

chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel

## [1] 7.814728

5. Keputusan Uji

Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0", "Tolak H0")
Keputusan

## [1] "Terima H0"

Interpretasi Pada taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi tiga arah antara jenis kelamin, tingkat pendidikan, dan pendapat mengenai kecepatana pelaayananan.

Catatan:

Model pengurang adalah model yang lebih lengkap (lebih banyak parameter, df lebih kecil), yaitu model saturated.
Derajat bebas dihitung dari selisih derajat bebas model homogenous dan saturated.
Keputusan berdasarkan perbandingan deviance model dengan chi-square tabel.

Rangkuman

Pengujian Ada Tidaknya Interaksi Tiga Arah (Saturated Model vs Homogenous Model)

Hipotesis
- \[H_0=\lambda_{ijk}^{XYZ}=0\] (Tidak ada interaksi tiga arah; model yang terbentuk adalah model homogenous)
- \[H_1=\lambda_{ijk}^{XYZ}\neq0\] (Ada interaksi tiga arah; model yang terbentuk adalah model saturated)
Tingkat Signifikansi
- α=5%
Statistik Uji
- ΔDeviance=Deviance model homogenous−Deviance model saturated
- =3.443845−3.347364e-10=3.443845
- db=db model homogenous−db model saturated=3-0 = 3
Daerah Penolakan
- Tolak H0 jika ΔDeviance>χ = 7.814
Keputusan
- Karena 3.443845 < 7.814 maka menerima H0
Interpretasi
- Pada taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi tiga arah antara jenis kelamin, pendidikan, dan pendapat mengenai kecepatanan pelayanan di RSHS

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.

\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{XZ} \]

# Conditional Association on X
model_conditional_X <- glm(counts ~ x.sex + y.fav + z.edu +
              x.sex*y.fav + x.sex*z.edu,
              family = poisson(link = "log"))
summary(model_conditional_X)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.edu + x.sex * y.fav + 
##     x.sex * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)   
## (Intercept)            -2.8332     1.3933  -2.034  0.04200 * 
## x.sexLaki               0.7538     1.8381   0.410  0.68175   
## y.favCepat              2.7726     1.0308   2.690  0.00715 **
## z.eduSMA                2.0794     1.0607   1.961  0.04994 * 
## z.eduSMP                1.3863     1.1180   1.240  0.21500   
## z.eduTinggi             1.3863     1.1180   1.240  0.21500   
## x.sexLaki:y.favCepat   -0.8267     1.2783  -0.647  0.51781   
## x.sexLaki:z.eduSMA     -0.1335     1.5059  -0.089  0.92934   
## x.sexLaki:z.eduSMP     -0.6931     1.6583  -0.418  0.67596   
## x.sexLaki:z.eduTinggi   0.4055     1.5546   0.261  0.79423   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 46.3151  on 15  degrees of freedom
## Residual deviance:  6.7886  on  6  degrees of freedom
## AIC: 55.35
## 
## Number of Fisher Scoring iterations: 6

Langkah-Langkah Pengujian

1. Hipotesis

\[H_0:\lambda_{jk}^{YZ}=0\](Tidak ada interaksi antara pendapat hukuman mati (Y) dan fundamentalisme (Z))
\[H_1:\lambda_{jk}^{YZ}\neq0\] (Ada interaksi antara pendapat hukuman mati (Y) dan fundamentalisme (Z))

2. Hitung Selisih Deviance

# Deviance of Model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance   # model_conditional_X: conditional on X, model_homogenous: homogenous
Deviance.model

## [1] 3.344759

3. Hitung Derajat Bebas

# Derajat bebas = db model homogenous - db model saturated 
derajat.bebas <- (model_conditional_X$df.residual - model_homogenous$df.residual) 
derajat.bebas

## [1] 3

4. Chi-Square Tabel

chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel

## [1] 7.814728

5. Keputusan Uji

Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan

## [1] "Terima"

Interpretasi : Pada Taraf Nyata 5%, belum cukup bukti untuk menolak H0 atau dengan kata lain tidak ada interaksi antara pendapat mengenai kecepatnaan pelayanan(Y) DAN Tingkat Pendidikan (Z) . Model yang terbentuk cukup hanya sampai 2 interaksi dengan X (Conditional onX), sehingga interaksi Y*Z tidak signifikan secara statistik.

UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON Y)

Model Conditional on Y

Model log-linear conditional pada Y memasukkan efek utama dan interaksi dua arah antara X dengan Y dan X dengan Z, tanpa interaksi antara X dengan Z maupun interaksi tiga arah.

\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{YZ} \]

# Conditional Association on Y
model_conditional_Y <- glm(counts ~ x.sex + y.fav + z.edu +
              x.sex*y.fav + y.fav*z.edu,
              family = poisson(link = "log"))
summary(model_conditional_Y)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.edu + x.sex * y.fav + 
##     y.fav * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)            -2.069e+01  1.088e+04  -0.002    0.998
## x.sexLaki               6.931e-01  1.225e+00   0.566    0.571
## y.favCepat              2.075e+01  1.088e+04   0.002    0.998
## z.eduSMA                2.028e+01  1.088e+04   0.002    0.999
## z.eduSMP                1.959e+01  1.088e+04   0.002    0.999
## z.eduTinggi             1.097e-07  1.539e+04   0.000    1.000
## x.sexLaki:y.favCepat   -8.267e-01  1.278e+00  -0.647    0.518
## y.favCepat:z.eduSMA    -1.841e+01  1.088e+04  -0.002    0.999
## y.favCepat:z.eduSMP    -1.867e+01  1.088e+04  -0.002    0.999
## y.favCepat:z.eduTinggi  1.609e+00  1.539e+04   0.000    1.000
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 46.3151  on 15  degrees of freedom
## Residual deviance:  4.9883  on  6  degrees of freedom
## AIC: 53.549
## 
## Number of Fisher Scoring iterations: 18

Langkah-Langkah Pengujian

1. Hipotesis

\[H_0:\lambda_{ik}^{XZ}=0\] (Tidak ada interaksi antara Jenis Kelamin (X) dan Tingkat Pedndidikan(Z))
\[H_1:\lambda_{ik}^{XZ}\neq0\] (Ada interaksi antara Jenis Kelamin (X) dan Tingkat Pedndidikan(Z))

2. Hitung Selisih Deviance

# Deviance of Model
Deviance.model <- model_conditional_Y$deviance - model_homogenous$deviance   # model_conditional_Y: conditional on Y, model_homogenous: homogenous
Deviance.model

## [1] 1.54446

3. Hitung Derajat Bebas

derajat.bebas <- (model_conditional_Y$df.residual - model_homogenous$df.residual) 
derajat.bebas

## [1] 3

4. Chi-Square Tabel

chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel

## [1] 7.814728

5. Keputusan Uji

Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan

## [1] "Terima"

Interpretasi : Pada Taraf Nyata 5%, belum cukup bukti untuk menolak H0 atau dengan kata lain tidak ada interaksi antara pendapat mengenai Jenis Kelamin (X) DAN Tingkat Pendidikan (Z) . Model yang terbentuk cukup hanya sampai 2 interaksi dengan Y (Conditional onY), sehingga interaksi X*Z tidak signifikan secara statistik.

UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON Z)

Model Conditional on Z

Model log-linear conditional pada Z memasukkan efek utama dan interaksi dua arah antara X dengan Y dan X dengan Z, tanpa interaksi antara X dengan Y maupun interaksi tiga arah

\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ik}^{XZ}+\lambda_{jk}^{YZ} \]

# Conditional Association on Z
model_conditional_Z <- glm(counts ~ x.sex + y.fav + z.edu +
              x.sex*z.edu + y.fav*z.edu,
              family = poisson(link = "log"))
summary(model_conditional_Z)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.edu + x.sex * z.edu + 
##     y.fav * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)            -2.030e+01  1.099e+04  -0.002    0.999
## x.sexLaki               2.011e-15  1.414e+00   0.000    1.000
## y.favCepat              2.030e+01  1.099e+04   0.002    0.999
## z.eduSMA                2.037e+01  1.099e+04   0.002    0.999
## z.eduSMP                1.990e+01  1.099e+04   0.002    0.999
## z.eduTinggi            -2.035e-01  1.547e+04   0.000    1.000
## x.sexLaki:z.eduSMA     -1.335e-01  1.506e+00  -0.089    0.929
## x.sexLaki:z.eduSMP     -6.931e-01  1.658e+00  -0.418    0.676
## x.sexLaki:z.eduTinggi   4.055e-01  1.555e+00   0.261    0.794
## y.favCepat:z.eduSMA    -1.843e+01  1.099e+04  -0.002    0.999
## y.favCepat:z.eduSMP    -1.869e+01  1.099e+04  -0.002    0.999
## y.favCepat:z.eduTinggi  1.590e+00  1.547e+04   0.000    1.000
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 46.3151  on 15  degrees of freedom
## Residual deviance:  4.3125  on  4  degrees of freedom
## AIC: 56.874
## 
## Number of Fisher Scoring iterations: 18

Langkah-Langkah Pengujian

1. Hipotesis

\[H_0:\lambda_{ij}^{XY}=0\] (Tidak ada interaksi antara Jenis Kelamin (X) dan Persepsi Kecepatan Pelayanan(Y))
\[H_1:\lambda_{ij}^{XY}\neq0\] (Ada interaksi antara Jenis Kelamin (X) dan Persepsi Kecepatan Pelayanan(Y))

2. Hitung Selisih Deviance

# Deviance of Model
Deviance.model <- model_conditional_Z$deviance - model_homogenous$deviance   # model_conditional_Z: conditional on Z, model_homogenous: homogenous
Deviance.model

## [1] 0.8686681

3. Hitung Derajat Bebas

derajat.bebas <- (model_conditional_Z$df.residual - model_homogenous$df.residual) 
derajat.bebas

## [1] 1

4. Chi-Square Tabel

chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel

## [1] 3.841459

5. Keputusan Uji

Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan

## [1] "Terima"

Interpretasi : Pada Taraf Nyata 5%, belum cukup bukti untuk menolak H0 atau dengan kata lain tidak ada interaksi antara pendapat mengenai Jenis Kelamin (X) DAN Tingkat Pendidikan (Z) . Model yang terbentuk cukup hanya sampai 2 interaksi dengan Z (Conditional onZ), sehingga interaksi X*Y tidak signifikan secara statistik.

Pemilihan Model Terbaik

Ringkasan Model Log Linear

Model	Parameter	Deviance	Jumlah Parameter	df	AIC
model_saturated	\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{XZ}+\lambda_{jk}^{YZ}+\lambda_{ijk}^{XYZ} \]	3.35e-10	16	0	60.561
model_homogenous	\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{XZ}+\lambda_{jk}^{YZ} \]	3.4438	13	3	58.005
model_conditional_X	\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{XZ} \]	6.7886	10	6	55.350
model_conditional_Y	\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{YZ} \]	4.9883	10	6	53.549
model_conditional_Z	\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ik}^{XZ}+\lambda_{jk}^{YZ} \]	4.3125	11	4	56.874

Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah

Interaksi	Pengujian	Δ deviance	Δ df	χ² Tabel	Keputusan	Keterangan
XYZ	Saturated vs Homogenous	3.4438	3	7.8147	Terima H0	tidak ada interaksi
YZ	Conditional on X vs Homogenous	3.3448	3	7.8147	Terima	tidak ada interaksi
XZ	Conditional on Y vs Homogenous	1.5445	3	7.8147	Terima	tidak ada interaksi
XY	Conditional on Z vs Homogenous	0.8687	1	3.8415	Terima	tidak ada interaksi

Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah

Dari hasil uji deviance, semua model tidak menolak H₀, yang berarti model-model yang lebih sederhana (tanpa interaksi tertentu) tidak berbeda signifikan dari model yang lebih kompleks. Maka kita akan gunakan alternatif pemilihan model terbaik berdasarkan AIC nya

Berdasarkan AIC didapatkan bahwa model condition on Y merupakan model terbaik karena memiliki nilai AIC terkecil yaitu 53.549. Sehingga model terbaik :

\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{YZ} \] Model terbaik adalah model log-linear tanpa interaksi tiga arah dan hanya memuat interaksi dua arah antara jenis kelamin dan persepsi terhadap kecepatan pelayanan di RSHS.

Model Terbaik

Model terbaik adalah model log-linear tanpa interaksi tiga arah dan hanya memuat interaksi dua arah antara jenis kelamin (X) dan persepsi terhadap kecepatan pelayanan di RSHS (Z).

\[ \lambda+\lambda_i^X+\lambda_j^Y+\lambda_k^Z+\lambda_{ij}^{XY}+\lambda_{ik}^{YZ} \]

# Model Terbaik
bestmodel <- glm(counts ~ x.sex + y.fav + z.edu +
                 x.sex*y.fav + y.fav*z.edu,
                 family = poisson(link = "log"))
summary(bestmodel)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.edu + x.sex * y.fav + 
##     y.fav * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)            -2.069e+01  1.088e+04  -0.002    0.998
## x.sexLaki               6.931e-01  1.225e+00   0.566    0.571
## y.favCepat              2.075e+01  1.088e+04   0.002    0.998
## z.eduSMA                2.028e+01  1.088e+04   0.002    0.999
## z.eduSMP                1.959e+01  1.088e+04   0.002    0.999
## z.eduTinggi             1.097e-07  1.539e+04   0.000    1.000
## x.sexLaki:y.favCepat   -8.267e-01  1.278e+00  -0.647    0.518
## y.favCepat:z.eduSMA    -1.841e+01  1.088e+04  -0.002    0.999
## y.favCepat:z.eduSMP    -1.867e+01  1.088e+04  -0.002    0.999
## y.favCepat:z.eduTinggi  1.609e+00  1.539e+04   0.000    1.000
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 46.3151  on 15  degrees of freedom
## Residual deviance:  4.9883  on  6  degrees of freedom
## AIC: 53.549
## 
## Number of Fisher Scoring iterations: 18

Interpretasi Koefisien Model Terbaik

# Interpretasi koefisien model terbaik
data.frame(
  koef = bestmodel$coefficients,
  exp_koef = exp(bestmodel$coefficients)
)

\(\exp(\lambda_{Intercept})\) = \(\exp(-20.68863) = 1.035239 \times 10^{-9}\) → nilai odds
Ini adalah odds dasar (baseline odds), tanpa memperhatikan jenis kelamin, preferensi, atau tingkat pendidikan.

\(\exp(\lambda_{x.sexLaki})\) = \(exp(0.6931472) = 2.000\) → nilai odds
Tanpa memperhatikan preferensi dan tingkat pendidikan, peluang seseorang berjenis kelamin laki-laki adalah 2 kali dibandingkan perempuan.
Atau, peluang perempuan adalah \(1/2 = 0.5\) kali dibandingkan laki-laki.

\(\exp(\lambda_{y.favCepat})\)= \(\exp(20.75317) = 1.030358 \times 10^9\) → nilai odds
Tanpa memperhatikan jenis kelamin dan tingkat pendidikan, peluang seseorang yang mendukung proses cepat adalah 1.03 miliar kali dibandingkan yang tidak mendukung.
\(\exp(\lambda_{z.eduSMA})\) = \(\exp(20.28317) = 6.439739 \times 10^8\) → nilai odds
Tanpa memperhatikan jenis kelamin dan preferensi, peluang seseorang dengan pendidikan SMA adalah 643 juta kali dibandingkan yang berpendidikan referensi (kategori referensi, misal: SD atau lainnya)

\(\exp(\lambda_{z.eduSMP})\) = \(\exp(19.59002) = 3.219870 \times 10^8\) → nilai odds
Tanpa memperhatikan jenis kelamin dan preferensi, peluang seseorang dengan pendidikan SMP adalah 321 juta kali dibandingkan kategori referensi.

\(\exp(\lambda_{z.eduTinggi})\) = \(\exp(1.096651 \times 10^{-7}) = 1.000\) → nilai odds
Tanpa memperhatikan jenis kelamin dan preferensi, peluang seseorang dengan pendidikan tinggi adalah 1 kali dibandingkan kategori referensi (tidak ada perbedaan yang berarti).

\(\exp(\lambda_{x.sexLaki:y.favCepat})\) = \(\exp(-0.8266786) = 0.4375\) → nilai odds ratio

Tanpa memperhatikan tingkat pendidikan, **odds mendukung proses cepat** jika dia laki-laki adalah **0.4375 kali** dibandingkan odds yang sama jika dia perempuan.\
Atau, perempuan memiliki peluang $1/0.4375 = 2.29$ kali lebih besar dibanding laki-laki.

\(\exp(\lambda_{y.favCepat:z.eduSMA})\) = \(\exp(-18.41137) = 1.009358 \times 10^{-8}\) → nilai odds ratio
Tanpa memperhatikan jenis kelamin, odds mendukung proses cepat pada pendidikan SMA adalah sangat kecil (sekitar \(1.009 \times 10^{-8}\) kali) dibandingkan kategori referensi.

\(\exp(\lambda_{y.favCepat:z.eduSMP})\) = \(\exp(-18.67373) = 7.764290 \times 10^{-9}\) → nilai odds ratio
Tanpa memperhatikan jenis kelamin, odds mendukung proses cepat pada pendidikan SMP adalah sangat kecil (sekitar \(7.76 \times 10^{-9}\) kali) dibandingkan kategori referensi.

\(\exp(\lambda_{y.favCepat:z.eduTinggi})\) = \(\exp(1.609438) = 5.000\) → nilai odds ratio
Tanpa memperhatikan jenis kelamin, odds mendukung proses cepat pada pendidikan Tinggi adalah 5 kali dibandingkan kategori referensi.

Nilai Dugaan Model Terbaik

# Fitted values dari model terbaik
data.frame(
Fund = z.edu,
sex = x.sex,
favor = y.fav,
counts = counts,
fitted = bestmodel$fitted.values
)

Perhitungan Manual Nilai Dugaan Model Terbaik

Secara manual, nilai ftted value diperoleh dengan cara sebagai berikut:

Perhitungan \(\hat{\mu}_{ijk}\)

Tinggi - Laki - Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0.6931472 + 20.75317 + 0 + 1.609438 + (-0.8266786)) = \exp(4.6462576) = 4.66667 \]

Tinggi - Laki - Tidak Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0.6931472 + 0 + 0 + 0 + 0) = \exp(-19.9954828) = 2.070477 \times 10^{-9} \]

Tinggi - Perempuan - Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0 + 20.75317 + 0 + 1.609438 + 0) = \exp(1.673978) = 5.33333 \]

Tinggi - Perempuan - Tidak Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0 + 0 + 0 + 0 + 0) = \exp(-20.68863) = 1.035239 \times 10^{-9} \]

SMA - Laki - Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0.6931472 + 20.75317 + 20.28317 + (-18.41137) + (-0.8266786)) = \exp(1.8038386) = 6.06667 \]

SMA - Laki - Tidak Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0.6931472 + 0 + 20.28317 + 0 + 0) = \exp(0.2876872) = 1.33333 \]

SMA - Perempuan - Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0 + 20.75317 + 20.28317 + (-18.41137) + 0) = \exp(1.1033366) = 6.93333 \]

SMA - Perempuan - Tidak Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0 + 0 + 20.28317 + 0 + 0) = \exp(-0.40546) = 0.66667 \]

SMP - Laki - Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0.6931472 + 20.75317 + 19.59002 + (-18.67373) + (-0.8266786)) = \exp(0.847982) = 2.33333 \]

SMP - Laki - Tidak Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0.6931472 + 0 + 19.59002 + 0 + 0) = \exp(-0.40546) = 0.66667 \]

SMP - Perempuan - Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0 + 20.75317 + 19.59002 + (-18.67373) + 0) = \exp(0.980752) = 2.66667 \]

SMP - Perempuan - Tidak Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0 + 0 + 19.59002 + 0 + 0) = \exp(-1.098617) = 0.33333 \]

SD - Laki - Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0.6931472 + 20.75317 + 0 + 0 + (-0.8266786)) = \exp(-0.069) = 0.93333 \]

SD - Laki - Tidak Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0.6931472 + 0 + 0 + 0 + 0) = \exp(-19.9954828) = 2.070477 \times 10^{-9} \]

SD - Perempuan - Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0 + 20.75317 + 0 + 0 + 0) = \exp(0.06454) = 1.06667 \]

SD - Perempuan - Tidak Cepat

\[ \hat{\mu} = \exp(-20.68863 + 0 + 0 + 0 + 0 + 0) = \exp(-20.68863) = 1.035239 \times 10^{-9} \]

Referensi

Fisher, R. A. (1935). The Design of Experiments. Oliver and Boyd.
Agresti, A. (2013). Categorical Data Analysis. John Wiley & Sons.
Mehta, C. R., & Patel, N. R. (1983). A Network Algorithm for Performing Fisher’s Exact Test in r × c Contingency Tables. Journal of the American Statistical Association, 78(382), 427-434.
Howell, D. C. (2012). Statistical Methods for Psychology (8th ed.). Cengage Learning
BPS Kabupaten Bandung. KABUPATEN BANDUNG DALAM ANGKA Bandung Regency in Figures 2024. vol. 41, Kabupaten Bandung, BPS Kabupaten Bandung, 2024,
BPS Kabupaten Bandung Barat. KABUPATEN BANDUNG BARAT DALAM ANGKA Bandung Barat Regency in Figures 2024. vol. 15, Kabupaten Bandung Barat, BPS Kabupaten Bandung Barat, 2024,
BPS Kota Bandung. KOTA BANDUNG DALAM ANGKA Bandung Municipality in Figures 2024. vol. 44, Bandung, BPS Kota Bandung, 2024,
BPS Kota Cimahi. KOTA CIMAHI DALAM ANGKA Cimahi Municipality in Figures 2024. vol. 2022, Kota Cimahi, BPS Kota Cimahi, 2024,
Aeberhard, S., and D. Coomans. “Wine.” UCI Machine Learning Repository, https://archive.ics.uci.edu/dataset/109/wine.