1 Pendahuluan

Dalam statistik, data kategori adalah jenis data yang dibagi ke dalam kelompok atau kategori dan bukan berupa angka yang bisa diukur. Data seperti ini sering digunakan di berbagai bidang untuk memahami pola, hubungan, dan tren yang tidak bisa dilihat langsung melalui angka. Contoh dari data kategori yaitu jenis kelamin (laki-laki/perempuan), status pernikahan (menikah/belum menikah), tingkat pendidikan (rendah/menengah/tinggi), serta preferensi konsumen (setuju/netral/tidak setuju). Analisis terhadap data kategori sangat penting untuk mengeksplorasi dan memahami informasi yang bersifat nominal atau ordinal.

Karena berbeda dari data numerik, data kategori perlu dianalisis dengan metode khusus, seperti tabel kontingensi, uji chi-square, regresi logistik, dan pendekatan berbasis probabilitas. Seiring berkembangnya teknologi, analisis data kategori juga makin sering digunakan dalam machine learning dan kecerdasan buatan.

1.1 Tujuan Analisis Data Kategori

Analisis Data Kategori memiliki banyak tujuan, diantaranya yaitu:

1.1.1 Melihat Frekuensi Data Tersebar dalam Tiap Kategori

Analisis data kategori memungkinkan peneliti untuk melihat bagaimana distribusi data tiap kategori. Misalkan terdapat responden dengan kategori usia muda, dewasa, dan tua. Peneliti dapat melihat sebaran usia dari responden yang ada.

1.1.2 Mengindentifikasi Pola atau Tren

Dalam analisis data kategori, identifikasi pola bertujuan untuk menemukan kecenderungan perbedaan antar kelompok dalam data. Misalnya, melalui visualisasi seperti diagram batang dan diagram lingkaran, kita dapat melihat kecenderungan responden memilih suatu merk dibandingkan dengan merk lain.

1.1.3 Mendukung dalam Pengambilan Keputusan

Dari tujuan yang telah dijelaskan sebelumnya, peneliti juga dimudahkan dalam pengambilan keputusan karena tren atau polanya sudah diketahui.

1.2 Definisi dan Ruang Lingkup Analisis Data Kategori

Analisis data kategorik merupakan suatu teknik statistik yang digunakan untuk menganalisis variabel yang bersifat kategorik, dimana data yang dianalisis bersifat nominal atau ordinal.

1.2.1 Definisi Nominal dan Ordinal

Data yang bersifat nominal tidak memiliki urutan, contohnya seperti jenis kelamin dan status menikah. Sedangkan data yang bersifat ordinal memiliki urutan, contohnya seperti tingkat pendidikan dan preferensi konsumen.

1.2.2 Data Biner dan Multikategori

Dalam analisis kategorik, data juga dapat dibedakan berdasarkan jumlah kategorinya. Data biner berarti data hanya memiliki dua kategori, misalnya ya/tidak. Sedangkan data multikategori berarti data tersebut memiliki lebih dari dua kategori.

1.3 Perbedaan Data Kategorik dan Data Numerik

Data Kategorik Data Kuantitatif
Bentuk Kategori Numerik
Sifat Kualitatif Kuantitatif
Skala Pengukuran Skala nominal/ordinal Skala Interval/ratio

1.4 Manfaat Analisis Data Kategorik dalam Berbagai Bidang

Analisis data kategori memiliki berbagai manfaat dalam berbagai bidang, diantaranya yaitu:

Bidang Tujuan Penerapan
Kesehatan Mengetahui hubungan antara kondisi pasien dan kategori risiko Melihat hubungan antara pola hidup (Sehat/Tidak Sehat) dengan penderita Diabetes
Sosial Menganalisis opini masyarakat Menganalisis kepuasan terhadap layanan publik (Sangat Puas – Tidak Puas)
Ekonomi Mengkategorikan jenis pekerjaan dan status ekonomi Melihat hubungan antara jenis pekerjaan dan tingkat pengeluaran

2 Metode dalam Analisis Data Kategori

Dalam analisis data kategori, metode yang digunakan bergantung kepada tujuan penelitian. Ada beberapa metode dalam analisis data kategori, diantaranya yaitu:

2.1 Uji Chi Square

Uji chi - square bertujuan untuk menguji hubungan antara dua variabel, contohnya yaitu menguji apakah ada hubungan antara jenis kelamin dan preferensi produk.

2.2 Uji Fisher’s Exact

Uji Fisher’s Exact bertujuan sama seperti uji chi - square, hanya saya uji ini digunakan ketika frekuensi data berjumlah kecil (<5).

2.3 Regresi Logistik

Regresi logistik pada analisis data kategorik bertujuan untuk memprediksi kemungkinan suatu kategori terjadi berdasarkan variabel prediktor. Misalnya memprediksi apakah seseorang siap menikah berdasarkan tingkat pendapatan, kesiapan mental dan fisik, serta kematangan emosional.

2.4 Tabulasi Silang (Cross tabulation)

Tabulasi silang atau cross tabulation memiliki tujuan untuk menampilkan distribusi frekuensi dari dua variabel dalam bentuk tabel. Misalnya menampilkan jumlah laki - laki dan perempuan dalam tiap kelompok penyakit.

2.5 Analisis Diskriminan (Untuk Data Kategorik Sebagai Target)

Analisis diskriminan bertujuan untuk mengklasifikasikan objek ke dalam kategori berdasarkan beberapa variabel, misalnya menentukan preferensi produk pelanggan jenis kelamin dan pendapatan.

Analisis data kategorik memiliki banyak metode untuk digunakan sesuai tujuan masing-masing. Oleh karena itu, peneliti harus memiliki tujuan dan pemahaman untuk bisa meneliti data kategorik menggunakan metode yang tepat.

Data

Tingkat pendidikan merupakan faktor yang berpengaruh terhadap partisipasi seseorang dalam angkatan kerja. Untuk memahami hubungan tersebut, analisis dilakukan menggunakan data dari Badan Pusat Statistik (BPS) tahun 2021–2023 yang mencakup jumlah penduduk usia kerja menurut pendidikan tertinggi dan status keikutsertaan dalam angkatan kerja. Data dikategorikan berdasarkan jenjang pendidikan (SD hingga Perguruan Tinggi) dan status (AK dan Bukan AK), kemudian dianalisis menggunakan tabel kontingensi, uji chi-square, ukuran asosiasi, dan uji Cochran-Mantel-Haenszel (CMH). Tujuan dari analisis ini adalah untuk menguji apakah terdapat hubungan signifikan antara pendidikan dan status kerja, serta apakah hubungan tersebut konsisten antar tahun.

##         Pendidikan   X2021   X2022 Total_2023 X2021.1 X2022.1 Persen_AK_2023
## 1         >= SD/MI 1461122 1284939    1144695   55.76   53.13          55.80
## 2          SMP/MTs 1563661 1465334    1532001   46.86   46.93          47.28
## 3       SMA/SMK/MA 3838202 3980426    3874632   66.95   67.08          68.78
## 4 Perguruan Tinggi 1403371 1596304    1771585   75.54   75.91          78.98
##        AK Bukan_AK
## 1  638740   505955
## 2  724330   807671
## 3 2664972  1209660
## 4 1399198   372387

3 Distribusi Probabilitas dalam Data Kategori

Variabel acak kategori adalah variabel yang hanya memiliki beberapa kategori diskrit sebagai hasilnya. Distribusi probabilitas dari variabel ini menggambarkan kemungkinan terjadinya setiap kategori.

3.1 Distribusi Bernoulli

Distribusi Bernoulli adalah distribusi probabilitas diskrit yang hanya memiliki dua kemungkinan hasil: sukses (1) atau gagal (0). Distribusi ini digunakan untuk percobaan yang hanya memiliki dua hasil, seperti lemparan koin (kepala atau ekor).

Distribusi Bernoulli memiliki satu parameter, yaitu:

  • \(p\): probabilitas sukses (nilai antara 0 dan 1)

Jika \(X\) adalah peubah acak yang mengikuti distribusi Bernoulli, maka:

\[ X \sim \text{Bernoulli}(p) \]

3.1.1 Rumus Distribusi Bernoulli

Probability Mass Function (PMF) dari distribusi Bernoulli didefinisikan sebagai:

\[ P(X = x) = \begin{cases} p & \text{jika } x = 1 \\\\ 1 - p & \text{jika } x = 0 \end{cases} \]

atau bisa ditulis dengan persamaan:

\[ P(X = x) = p^x (1 - p)^{1 - x}, \quad x \in \{0, 1\} \]

3.1.2 Perhitungan dengan R

library(knitr)
library(kableExtra)
library(dplyr)
# Simulasi Bernoulli (masuk angkatan kerja = 1, tidak masuk = 0)
set.seed(123)  
bernoulli_SD <- rbinom(n = 10, size = 1, prob = 0.558)
bernoulli_SMP <- rbinom(n = 10, size = 1, prob = 0.473)
bernoulli_SMA <- rbinom(n = 10, size = 1, prob = 0.688)
bernoulli_PT <- rbinom(n = 10, size = 1, prob = 0.7898)
# hasil
bernoulli_SD
##  [1] 1 0 1 0 0 1 1 0 1 1
bernoulli_SMP
##  [1] 1 0 1 1 0 1 0 0 0 1
bernoulli_SMA
##  [1] 0 0 1 0 1 0 1 1 1 1
bernoulli_PT
##  [1] 0 0 1 0 1 1 1 1 1 1

Interpretasi: Hasil simulasi distribusi Bernoulli untuk masing-masing jenjang pendidikan menunjukkan variasi dalam peluang keberhasilan (ditandai dengan nilai 1) dan kegagalan (nilai 0). Pada kelompok SD, terdapat 6 keberhasilan dari 10 percobaan, yang mengindikasikan tingkat keberhasilan sebesar 60%. Untuk kelompok SMP, keberhasilan terjadi sebanyak 6 kali juga, menunjukkan proporsi yang sama. Pada kelompok SMA, terdapat peningkatan keberhasilan menjadi 7 dari 10 percobaan (70%), sedangkan kelompok perguruan tinggi (PT) menunjukkan tingkat keberhasilan tertinggi, yaitu 8 dari 10 percobaan atau 80%. Secara umum, tren ini menggambarkan bahwa semakin tinggi jenjang pendidikan, semakin besar pula kemungkinan keberhasilan dalam konteks percobaan biner yang disimulasikan dengan distribusi Bernoulli.

3.2 Distribusi Binomial

Distribusi Binomial adalah distribusi probabilitas diskrit yang menggambarkan jumlah keberhasilan dalam sejumlah percobaan Bernoulli yang independen dan identik. Setiap percobaan memiliki dua hasil: sukses atau gagal, dan probabilitas sukses tetap di setiap percobaan.

Jika \(X\) adalah peubah acak yang mengikuti distribusi binomial, maka:

\[ X \sim \text{Binomial}(n, p) \]

dengan:

  • \(n\): jumlah percobaan
  • \(p\): probabilitas sukses pada setiap percobaan

3.2.1 Rumus Distribusi Binomial

Probability Mass Functiion (PMF) dari distribusi binomial yaitu:

\[ P(X = k) = \binom{n}{k} p^k (1 - p)^{n - k}, \quad k = 0, 1, 2, ..., n \]

Dengan:

  • \(\binom{n}{k}\) adalah koefisien binomial = \(\dfrac{n!}{k!(n-k)!}\)
  • \(k\): jumlah keberhasilan

3.2.2 Perhitungan di R

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

3.3 Distribusi Multinomial

Distribusi Multinomial adalah perluasan dari distribusi Binomial yang mamiliki lebih dari dua hasil. Distribusi ini digunakan untuk menghitung probabilitas hasil dari sejumlah percobaan independen dimana tiap percobaan memiliki lebih dari dua kemungkinan hasil.

Jika terdapat \(k\) kategori dan setiap percobaan memiliki probabilitas:

\[ p_1, p_2, ..., p_k \quad \text{dengan} \quad \sum_{i=1}^k p_i = 1 \]

dan dilakukan \(n\) percobaan, maka jumlah hasil dari setiap kategori mengikuti distribusi multinomial sebagai berikut:

\[ (X_1, X_2, ..., X_k) \sim \text{Multinomial}(n; p_1, p_2, ..., p_k) \]

3.3.1 Rumus Distribusi Multinomial

Probability Mass Functiion dari distribusi Multinomial yaitu:

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

Dengan syarat:

  • \(\sum_{i=1}^k x_i = n\)
  • \(\sum_{i=1}^k p_i = 1\)

3.3.2 Perhitungan dengan R

set.seed(123)
multinomial <- rmultinom(n = 1, size = 15, prob = c(0.2, 0.8, 0.6))
multinomial
##      [,1]
## [1,]    1
## [2,]    7
## [3,]    7

3.4 Distribusi Poisson

Distribusi Poisson adalah distribusi probabilitas diskrit yang menggambarkan jumlah kejadian dalam suatu interval waktu tertentu, dengan asumsi bahwa:

  • Kejadian terjadi secara acak,
  • Tidak saling bergantung satu sama lain,
  • Rata-rata kejadian per interval konstan.

Jika \(X\) adalah peubah acak yang mengikuti distribusi Poisson dengan parameter \(\lambda\), maka:

\[ X \sim \text{Poisson}(\lambda) \]

dimana:

  • \(\lambda\): rata-rata jumlah kejadian dalam suatu interval.

3.4.1 Rumus Distribusi Poisson

Probability Mass Functiion (PMF) dari distribusi Poisson adalah:

\[ P(X = k) = \frac{e^{-\lambda} \lambda^k}{k!}, \quad k = 0, 1, 2, ... \]

Sifat-sifat penting:

  • Mean: \(\mu = \lambda\)
  • Varians: \(\sigma^2 = \lambda\)

3.4.2 Perhitungan dengan R

set.seed(123)
poisson <- rpois(15, lambda = 5) 
poisson
##  [1] 4 7 4 8 9 2 5 8 5 5 9 5 6 5 2

4 Desain Sampling dalam Analisis Data Kategori

Desain sampling adalah tahap awal dalam proses penelitian, terkhusus dalam pengumpulan data. Tujuannya adalah agar data yang dikumpulkan dapat digunakan untuk menarik kesimpulan yang akurat tentang populasi.

Dalam analisis data kategori, desain sampling merujuk pada metode yang digunakan untuk memilih sampel dari populasi agar data yang diperoleh dapat mewakili keseluruhan populasi dengan baik. Desain sampling sangatlah penting karena pemilihan sampel akan memengaruhi validitas hasil penelitian.

Desain sampling dalam analisis data kategori dapat diklasifikasikan menjadi dua pendekatan utama, yaitu prospective sampling dan retrosprective sampling.

4.1 Prospective Sampling

Prospective sampling adalah metode pengambilan sampel dimana peneliti mengikuti partisipan dari waktu sekarang ke masa depan untuk mengamati kejadian atau hasil tertentu. Prospective sampling biasanya digunakan dalam studi kohort prospektif. Beberapa jenis desain sampling dalam metode ini meliputi:

4.1.1 Eksperimen

Dalam studi eksperimen, peneliti menentukan intervensi, yaitu peserta secara acak dibagi ke dalam kelompok perlakuan (menerima intervensi) dan kelompok kontrol (tidak menerima intervensi atau menerima plasebo). Teknik sampling yang digunakan meliputi:

  • Simple Random Sampling (SRS) : Setiap individu memiliki probabilitas yang sama untuk dipilih.
  • Stratified Random Sampling : Populasi dibagi menjadi strata berdasarkan karakteristik tertentu, lalu sampel diambil secara acak dari tiap strata.
  • Clustered Sampling : Populasi dibagi menjadi cluster, kemudian beberapa cluster dipilih secara acak.

4.1.2 Studi Kohort

Studi kohort adalah jenis penelitian dimana peneliti mengamati peserta tanpa intervensi langsung dimana peserta dikelompokkan berdasarkan paparan alami terhadap faktor risiko (misalnya, merokok atau tidak merokok. Jenis sampling yang umum dalam studi yaitu:

  • Systematic Sampling : Subjek dipilih berdasarkan interval tertentu dari populasi
  • Census Sampling : Seluruh anggota dalam populasi tertentu diikutsertakan dalam penelitian

4.2 Retrospective Sampling

Retrospective sampling adalah metode pengambilan sampel dalam penelitian dimana peneliti melihat ke belakang untuk mengidentifikasi subjek berdasarkan outcome yang sudah terjadi, lalu menelusuri paparan atau faktor penyebabnya sebelumnya.

4.2.1 Studi Kasus-Kontrol

Case-Control Study adalah salah satu jenis studi observasional yang digunakan untuk menyelidiki hubungan antara faktor risiko dan suatu outcome (misalnya, penyakit). Berbeda dengan studi kohort, studi ini dimulai dari outcome dan bekerja mundur untuk menentukan paparan sebelumnya. Teknik sampling yang sering digunakan meliputi:

  • Purposive Sampling : Pemilihan sampel berdasarkan karakteristik yang relevan dengan tujuan penelitian.
  • Snowball Sampling : Responden awal membantu merekrut subjek lain yang memiliki karakteristik serupa.

4.2.2 Studi Kohort Retrospektif

Dalam studi kohort retrospektif, data historis digunakan untuk mengelompokkan individu berdasarkan paparan kemudian peneliti menganalisis hasil yang terjadi. Teknik sampling yang sering digunakan meliputi :

  • Quota Sampling : Sampel dipilih untuk mencerminkan proporsi tertentu yang ingin diteliti dalam populasi.
  • Convenience Sampling : Subjek dipilih berdasarkan ketersediaan data yang sudah ada.

5 Tabel Kontingensi 2 x 2

Tabel kontingensi 2x2 adalah tabel yang digunakan untuk menyajikan data kategorik dari dua variabel yang masing-masing memiliki dua kategori. Tabel ini berguna untuk melihat hubungan atau asosiasi antara dua variabel, misalnya perhitungan risiko relatif (relative risk), perbedaan risiko (risk difference) dan odds ratio.

Struktur umum Tabel Kontingensi 2x2:

Kategori B1 Kategori B2 Total
Kategori A1 a b a + b
Kategori A2 c d c + d
Total a + c b + d n

5.1 Distribusi Peluang dalam Tabel Kontingensi 2x2

5.1.1 Peluang Bersama

Peluang bersama adalah probabilitas kedua variabel terjadi secara bersamaan dalam suatu sel tabel kontingensi:

\[ P(A_i, B_j) = \frac{n_{ij}}{n} \]

5.1.2 Peluang Marginal

Peluang marginal adalah probabilitas suatu variabel tanpa mempertimbangkan variabel lainnya.

  • Peluang marginal baris:

\[ P(A_i) = \frac{n_{i.}}{n} \]

  • Peluang marginal kolom:

\[ P(B_j) = \frac{n_{.j}}{n} \]

5.1.3 Peluang Bersyarat

Peluang bersyarat adalah probabilitas suatu kejadian terjadi dengan syarat kejadian lain telah terjadi:

\[ P(B_j | A_i) = \frac{P(A_i, B_j)}{P(A_i)} = \frac{n_{ij}}{n_{i.}} \]

5.2 Ukuran Asosiasi dalam Data Kategori 2 x 2

Ukuran asosiasi dalam analisis data kategorik digunakan untuk mengukur kekuatan hubungan antara dua variabel yang bersifat kategorik. Ukuran asosiasi ini penting karena digunakan untuk mengetahui apakah hubungan yang ada signifikan dan seberapa kuat hubungan tersebut. Berikut adalah beberapa ukuran asosiasi:

5.2.1 Risk Difference (RD)

Risk Difference mengukur selisih risiko antara dua kelompok. Nilai RD menunjukkan perbedaan proporsi kejadian antara kelompok yang terpapar dan tidak terpapar.

Rumus:

\[ RD = \frac{a}{a + b} - \frac{c}{c + d} \] Dengan

  • Jika 𝑅𝐷>0, maka risiko kejadian lebih tinggi di Kelompok 1 dibandingkan Kelompok 2

  • Jika 𝑅𝐷<0, maka risiko kejadian lebih rendah di Kelompok 1 dibandingkan Kelompok 2

  • Jika 𝑅𝐷=0, maka tidak ada perbedaan risiko antara dua kelompok.

5.2.2 Relative Risk (RR)

Relative Risk membandingkan kemungkinan terjadinya suatu kejadian pada kelompok yang terpapar terhadap kelompok yang tidak terpapar.

Rumus:

\[ RR = \frac{a / (a + b)}{c / (c + d)} \] Dengan

  • Jika 𝑅𝑅>1, maka kejadian lebih sering terjadi di Kelompok 1 dibandingkan Kelompok 2.

  • Jika 𝑅𝑅<1, maka kejadian lebih jarang terjadi di Kelompok 1 dibandingkan Kelompok 2.

  • Jika 𝑅𝑅=1, maka tidak ada perbedaan risiko antara dua kelompok.

5.2.3 Odds Ratio (OR)

Odds Ratio digunakan untuk membandingkan peluang terjadinya suatu kejadian antara dua kelompok.

Rumus:

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

Dengan

  • Jika 𝑂𝑅 >1, maka peluang kejadian lebih besar di Kelompok 1 dibandingkan Kelompok 2.

  • Jika 𝑂𝑅 <1, maka peluang kejadian lebih kecil di Kelompok 1 dibandingkan Kelompok 2.

  • Jika 𝑂𝑅 =1, maka tidak ada perbedaan peluang kejadian antara dua kelompok.

5.2.4 Contoh Penerapan dalam Studi Kasus

5.2.4.1 Contoh Perhitungan Manual

  • Risk Difference \[ RD = \frac{a}{a + b} - \frac{c}{c + d} = \frac{1601}{164128} - \frac{510}{412878} = 0.0098 - 0.0012 = 0.0086 \]

  • Relative Risk \[ RR = \frac{a / (a + b)}{c / (c + d)} = \frac{1601/164128}{510/412878} = 0.0098 / 0.0012 = 8.17 \]

  • Odds Ratio \[ OR = \frac{a \cdot d}{b \cdot c} = \frac{1601 \cdot 412368}{162527 \cdot 510} = \frac{660201168}{82988770} = 7.96 \]

5.2.4.2 Contoh Perhitungan dengan R
## [1] 0.1020622

Interpretasi : Nilai Risk Difference (RD) sebesar 0.1021 menunjukkan bahwa terdapat selisih risiko sebesar 10,21% antara dua kelompok yang dibandingkan (misalnya, antara lulusan Perguruan Tinggi dan SMA/SMK/MA) dalam hal menjadi angkatan kerja (AK). Dengan kata lain, kelompok yang berada di tingkat pendidikan lebih tinggi memiliki peluang menjadi angkatan kerja sekitar 10% lebih besar dibanding kelompok dengan pendidikan yang lebih rendah. Nilai RD yang positif ini mengindikasikan adanya keunggulan relatif pada kelompok pertama dalam kaitannya dengan status pekerjaan.

  • RD = 0.0086 → Perbedaan risiko kecelakaan fatal sebesar 0.86% lebih tinggi pada pengendara yang tidak pakai sabuk.

  • RR = 8.17 → Pengendara yang tidak memakai sabuk memiliki 8.17 kali risiko lebih besar dalam mengalami kecelakaan fatal.

  • OR = 7.96 → Peluang kecelakaan fatal 7.96 kali lebih besar pada pengendara yang tidak memakai sabuk.

6 Inferensi Tabel Kontingensi Dua Arah

Inferensi pada statistik merujuk pada proses penarikan kesimpulan tentang suatu populasi berdasarkan data dari sampel. Dalam konteks tabel kontingensi dua arah, inferensi digunakan untuk mengevaluasi dan memahami hubungan antara dua variabel kategorikal yang ditampilkan dalam bentuk tabel matriks.

Tabel kontingensi sendiri menyajikan distribusi frekuensi dari kombinasi kategori dua variabel. Tujuan utama dari inferensi pada tabel kontingensi adalah untuk menganalisis keterkaitan antarvariabel dengan melakukan proses estimasi dan pengujian hipotesis. Secara umum, proses inferensi dalam tabel kontingensi dua arah mencakup dua pendekatan utama yaitu Estimasi dan Pengujian

6.1 Estimasi

Estimasi berfungsi untuk memperkirakan parameter populasi berdasarkan sampel. Estimasi dibagi menjadi dua, yaitu:

6.1.1 Estimasi Titik

Estimasi titik digunakan untuk menentukan satu nilai spesifik sebagai perkiraan terbaik dari parameter populasi \[ \hat{p} = \frac{x}{n} \] Dengan:

  • \(\hat{p}\): Estimasi titik proporsi
  • \(x\): Jumlah individu dalam kategori tertentu
  • \(n\): Total jumlah individu dalam sampel

6.1.2 Estimasi Interval

Estimasi interval berfungsi untuk memberikan rentang nilai yang diyakini mengandung parameter populasi dengan tingkat kepercayaan tertentu

\[ \hat{p} \pm Z_{\alpha/2} \sqrt{ \hat{p}(1 - \hat{p}) } \] Dengan:

  • \(Z_{\alpha/2}\): Nilai dari distribusi normal standar untuk tingkat kepercayaan tertentu
  • \(\hat{p}\): Estimasi titik proporsi
  • \(n\): Ukuran sampel

6.2 Uji Hipotesis

6.2.1 Uji Proporsi

Uji ini digunakan untuk membandingkan proporsi kejadian antara dua kelompok dalam tabel kontingensi, terutama untuk menentukan apakah terdapat perbedaan yang signifikan dalam proporsi kejadian antara dua kelompok yang berbeda. Tabel kontingensi 2×2 memiliki struktur sebagai berikut:

Kejadian (+) Tidak Kejadian (-) Total (n)
Group 1 n11 n12 n1.
Group 2 n21 n22 n2.
Total n.1 n.2 n

Untuk menguji hipotesis bahwa tidak ada perbedaan proporsi antara dua kelompok, kita menggunakan uji z dua proporsi, dengan hipotesis:

- Hipotesis Nol (H0) : tidak ada perbedaan proporsi antara dua kelompok

- Hipotesis Alternatif (H1) : terdapat perbedaan proporsi antara dua kelompok

Estimasi proporsi dalam masing - masing kelompok diberikan oleh :

\[ \hat{p}_1 = \frac{a}{a + b} \]

\[ \hat{p}_2 = \frac{c}{c + d} \]

Estimasi proporsi gabungan

\[ \hat{p} = \frac{x_1 + x_2}{n_1 + n_2} \]

Statistik uji untuk uji proporsi dua sampel :

\[ Z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{\hat{p}(1 - \hat{p})\left(\frac{1}{n_1} + \frac{1}{n_2}\right)}} \]

Setelah nilai Z diketahui, kita bisa melihat keputusan apakah tolak H0 atau terima H0 dengan cara apabila |Z| lebih besar dari nilai kritis tertentu untuk tingkat signifikansi tertentu, maka H0 ditolak.

6.2.1.2. Perhitungan dengan R
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  x out of n
## X-squared = 62306, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.100859 0.102373
## sample estimates:
##    prop 1    prop 2 
## 0.7900112 0.6883952

Interpretasi : Jika kita menggunakan taraf signifikansi sebesar 5%, kita menolak H0 karena nilai |Z| yaitu 264.49 > 1.96, yang berarti bahwa ada perbedaan proporsi yang signifikan antara kelompok lulusan Perguruan Tinggi dan SMA/SMK/MA dalam hal status angkatan kerja.

6.2.2 Uji Asosiasi

Uji asosiasi dalam tabel kontingensi 2×2 bertujuan untuk mengukur hubungan antara dua variabel kategori. Tiga ukuran utama dalam uji asosiasi adalah: 1. Risk Diference (RD) – Mengukur selisih risiko absolut antara dua kelompok. 2. Relative Risk (RR) – Mengukur perbandingan risiko antara dua kelompok. 3. Odds Ratio (OR) – Mengukur perbandingan odds antara dua kelompok. Hipotesis Uji dalam Tabel Kontingensi 2 × 2 Untuk setiap uji asosiasi, hipotesis yang diuji adalah: • Hipotesis Nol (H0): Tidak ada asosiasi antara dua variabel. • Hipotesis Alternatif (H1): Terdapat asosiasi antara dua variabel.

6.2.2.1 Risk Difference (RD)

Risk Difference mengukur perbedaan absolut dalam probabilitas kejadian antara dua kelompok:

\[ RD = \left( \frac{n_{11}}{n_{1.}} \right) - \left( \frac{n_{21}}{n_{2.}} \right) \]

Standard Error untuk RD:

\[ SE(RD) = \sqrt{ \frac{\hat{p}_1 (1 - \hat{p}_1)}{n_{1.}} + \frac{\hat{p}_2 (1 - \hat{p}_2)}{n_{2.}} } \]

Statistik Uji Z:

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

6.2.2.2 Relative Risk (RR)

Relative Risk membandingkan kemungkinan kejadian antara dua kelompok:

\[ RR = \frac{n_{11} / n_{1.}}{n_{21} / n_{2.}} \]

Standard Error untuk log(RR):

\[ SE(\ln RR) = \sqrt{ \frac{1}{n_{11}} - \frac{1}{n_{1.}} + \frac{1}{n_{21}} - \frac{1}{n_{2.}} } \]

Statistik Uji Z:

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

Odds Ratio (OR)

Odds Ratio membandingkan odds kejadian antara dua kelompok:

\[ OR = \frac{n_{11} \times n_{22}}{n_{12} \times n_{21}} \]

Standard Error untuk log(OR):

\[ SE(\ln OR) = \sqrt{ \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} } \]

Statistik Uji Z:

\[ Z_{OR} = \frac{\ln OR}{SE(\ln OR)} \] Perhitungan menggunakan R

## $Risk_Difference
## [1] 0.1020622
## 
## $SE_RD
## [1] 0.0003858884
## 
## $Z_RD
## [1] 264.4863
## 
## $Relative_Risk
## [1] 1.148261
## 
## $SE_Ln_RR
## [1] 0.0005162791
## 
## $Z_RR
## [1] 267.779
## 
## $Odds_Ratio
## [1] 1.707546
## 
## $SE_Ln_OR
## [1] 0.002147789
## 
## $Z_OR
## [1] 249.12

Interpretasi: Berdasarkan hasil perhitungan, dapat disimpulkan bahwa terdapat hubungan yang signifikan antara jenjang pendidikan dan status angkatan kerja, khususnya saat membandingkan lulusan Perguruan Tinggi dengan lulusan SMA/SMK/MA:

  • Risk Difference (RD) sebesar 0.1021 menunjukkan bahwa proporsi lulusan Perguruan Tinggi yang masuk angkatan kerja lebih tinggi sekitar 10.2% dibanding lulusan SMA/SMK/MA.
  • Relative Risk (RR) sebesar 1.1483 mengindikasikan bahwa lulusan Perguruan Tinggi memiliki kemungkinan 14.8% lebih besar untuk masuk angkatan kerja dibanding lulusan SMA/SMK/MA.
  • Odds Ratio (OR) sebesar 1.7075 menunjukkan bahwa peluang lulusan Perguruan Tinggi untuk bekerja 1.7 kali lebih besar dibandingkan lulusan SMA/SMK/MA.
  • Nilai Z untuk RD, RR, dan OR semuanya sangat besar (di atas 249), dan didasarkan pada nilai Standard Error yang kecil, ini mengarah pada kesimpulan bahwa hasil tersebut sangat signifikan secara statistik.

Dengan demikian, perbedaan tingkat keterlibatan dalam angkatan kerja antara lulusan Perguruan Tinggi dan SMA/SMK/MA tidak hanya terlihat secara praktis, tetapi juga kuat secara statistik.

6.2.3. Uji Independensi

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

6.2.3.1. Uji Chi Square

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

Rumus Chi-Square

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

di mana:

  • \(O\) adalah nilai observasi dalam tabel kontingensi.
  • \(E\) adalah nilai yang diharapkan, dihitung sebagai:

\[ E_{ij} = \frac{R_i \times C_j}{N} \] dengan:

  • \(R_i\) = total baris ke-\(i\),
  • \(C_j\) = total kolom ke-\(j\),
  • \(N\) = total sampel.

Perhitungan menggunakan R

##                       AK Bukan_AK
## >= SD/MI          638740   505955
## SMP/MTs           724330   807671
## SMA/SMK/MA       2664972  1209660
## Perguruan Tinggi 1399198   372387
## 
##  Pearson's Chi-squared test
## 
## data:  data_chisq
## X-squared = 431600, df = 3, p-value < 2.2e-16

Interpretasi: Dengan taraf signifikansi 5%, didapatkan p-value < 0.05, maka kita menolak H0 (hipotesis nol). Artinya, terdapat hubungan yang signifikan secara statistik antara tingkat pendidikan (Perguruan Tinggi vs. SMA/SMK/MA) dan status angkatan kerja, di mana proporsi lulusan Perguruan Tinggi yang bekerja (0.7900) lebih tinggi secara signifikan dibandingkan lulusan SMA/SMK/MA (0.6884).

6.2.3.2. Uji Partisi Chi Square

Partisi Chi-Square Partisi Chi-Square digunakan untuk mengidentifikasi kategori mana dalam tabel kontingensi yang bertanggung jawab atas hubungan yang signifikan. Jika uji Chi-Square pada tabel kontingensi I×Jsignifikan, maka partisi Chi-Square memungkinkan kita untuk menguraikan efek hubungan dalam subkelompok yang lebih kecil.

6.2.3.2.1. Contoh Perhitungan

Perhitungan Dengan R:

##                        AK Bukan AK
## Pendidikan Rendah 1363070  1313626
## Pendidikan Tinggi 4064170  1582047
## $Chi_Square
## X-squared 
##  354884.7 
## 
## $P_Value
## [1] 0
## 
## $Decision
## [1] "Tolak H0: Ada hubungan"

Interpretasi: Karena nilai \(\chi^2 = 354{,}884.7\) sangat besar dan melebihi nilai kritis untuk \(df = 3\) pada taraf signifikansi 0.05, serta p-value = 0 (lebih kecil dari 0.05), maka kita menolak \(H_0\). Artinya, terdapat hubungan yang signifikan secara statistik antara tingkat pendidikan dan status angkatan kerja.

6.2.3.3. Uji Likelihood Ratio

Uji Likelihood Ratio (G²) adalah alternatif dari uji chi-square yang digunakan untuk menguji hipotesis independensi dalam tabel kontingensi \(I \times J\). Statistik uji ini diberikan oleh:

\[ G^2 = 2 \sum_{i} \sum_{j} n_{ij} \ln \left( \frac{n_{ij}}{\hat{\mu}_{ij}} \right) \]

Dimana:

- \(n_{ij}\) adalah frekuensi observasi dalam tabel kontingensi.

  • \(\hat{\mu}_{ij} = n \cdot p_i \cdot p_j\) adalah frekuensi yang diharapkan.

  • \(p_i\) adalah proporsi untuk baris \(i\).

- \(p_j\) adalah proporsi untuk kolom \(j\).

  • \(n\) adalah total frekuensi.
6.2.3.3.1. Contoh Perhitungan

Perhitungan Dengan R

## $G2
## [1] 431832.5
## 
## $df
## [1] 3
## 
## $Critical_Value
## [1] 7.814728
## 
## $Decision
## [1] "Tolak H0: Ada hubungan"

Interpretasi: Karena \(G^2 = 431{,}832.5 > 7.814728\), maka kita menolak hipotesis nol. Artinya, terdapat hubungan yang signifikan secara statistik antara tingkat pendidikan dan status angkatan kerja.

# Pastikan package-nya sudah terpasang
library(DescTools)
# Data dari data_ak
# Baris:
# 1 = SMA
# 2 = PT
# Kolom:
# 1 = AK
# 2 = Bukan AK
# Layer (strata): SD dan SMP (asumsi)
stratum1 <- matrix(c(data_ak$AK[3], data_ak$Bukan_AK[3], data_ak$AK[4], data_ak$Bukan_AK[4]), nrow=2, byrow=TRUE)

# Misal layer 2 bisa kita duplikat contoh aja (jika kamu hanya punya satu strata bisa buat dummy strata untuk contoh)
stratum2 <- matrix(c(3874632, 1800000, 1771585, 470000), nrow=2, byrow=TRUE)

# Gabungkan jadi array
pendidikan_array <- array(c(stratum1, stratum2),
  dim = c(2, 2, 2),
  dimnames = list(
    Pendidikan = c("SMA", "PT"),
    Status = c("AK", "Bukan AK"),
    Strata = c("Strata 1", "Strata 2")
  )
)

# Tampilkan arraynya
pendidikan_array
## , , Strata = Strata 1
## 
##           Status
## Pendidikan      AK Bukan AK
##        SMA 2664972  1209660
##        PT  1399198   372387
## 
## , , Strata = Strata 2
## 
##           Status
## Pendidikan      AK Bukan AK
##        SMA 3874632  1800000
##        PT  1771585   470000
# Jalankan uji
breslow_result <- BreslowDayTest(pendidikan_array)

# Tampilkan hasil
print(breslow_result)
## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  pendidikan_array
## X-squared = 85.747, df = 1, p-value < 2.2e-16

Interpretasi: Hasil Uji Breslow-Day pada array pendidikan_array menunjukkan nilai X-squared sebesar 85.747 dengan derajat bebas (df) = 1 dan p-value < 2.2e-16. Karena p-value jauh lebih kecil dari 0.05, maka kita menolak hipotesis nol (H₀) yang menyatakan bahwa odds ratio homogen di seluruh strata. Artinya, terdapat perbedaan yang signifikan pada odds ratio antara strata 1 dan strata 2, atau dengan kata lain, hubungan antara tingkat pendidikan dan status AK tidak konsisten di antara kedua strata tersebut. Ini mengindikasikan adanya interaksi antara faktor pendidikan dan strata terhadap peluang bekerja di sektor AK.

7. Tabel Kontingensi Tiga Arah

Bab ini membahas analisis tabel kontingensi tiga arah, mencakup tabel parsial dan marginal, ukuran asosiasi, uji independensi bersyarat, dan Simpson’s Paradox. Fokus utamanya adalah menganalisis hubungan antara dua variabel kategori dengan mempertimbangkan pengaruh variabel ketiga (kovariat). Tabel ini penting dalam banyak konteks, seperti epidemiologi, sosial, atau hukum.

Tabel tiga arah dapat dikembangkan menjadi:

Tabel Parsial: Mengisolasi efek kovariat (Z), untuk melihat hubungan bersyarat antara X dan Y.

Tabel Marginal: Mengabaikan Z, memberikan gambaran umum namun bisa menyesatkan (misalnya, dalam Simpson’s Paradox).

Penggunaan tabel marginal cocok untuk melihat asosiasi secara keseluruhan, tetapi bisa menyembunyikan informasi penting. Oleh karena itu, analisis tabel parsial sering lebih disarankan untuk hasil yang akurat dalam penelitian ilmiah. Jadikan Tahun untuk strata (dibagi jadi 3 strata)

7.1 Tabel Parsial dan Marginal

Tabel parsial adalah tabel yang mengelompokkan variabel \(X\) dan \(Y\) berdasarkan setiap level dari variabel ketiga \(Z\),
sedangkan tabel marginal adalah tabel yang mengabaikan \(Z\), dengan menjumlahkan data dari semua level \(Z\).

Struktur Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah dapat digambarkan secara umum sebagai berikut:

Z = 1 Z = 2 … Z = k
Y = 1 Y = 2 … Y = J Total Y = 1 Y = 2 … Y = J Total … Y = 1 … Y = J
X = 1 \(n_{111}\) \(n_{121}\) … \(n_{1+1}\) \(n_{112}\) \(n_{122}\) … \(n_{1+2}\) … \(n_{11k}\) … \(n_{1+k}\)
X = 2 \(n_{211}\) \(n_{221}\) … \(n_{2+1}\) \(n_{212}\) \(n_{222}\) … \(n_{2+2}\) … \(n_{21k}\) … \(n_{2+k}\)
… … … … … … … … … … … … …
X = I \(n_{I11}\) \(n_{I21}\) … \(n_{I+1}\) \(n_{I12}\) \(n_{I22}\) … \(n_{I+2}\) … \(n_{I1k}\) … \(n_{I+k}\)
Total \(n_{+11}\) \(n_{+21}\) … \(n_{++1}\) \(n_{+12}\) \(n_{+22}\) … \(n_{++2}\) … \(n_{+1k}\) … \(n_{++k}\)

Perhitungan Menggunakan R

##   Tahun       Pendidikan   Total      AK Bukan_AK
## 1  2021       SMA/SMK/MA 3838202 2569676  1268526
## 2  2021 Perguruan Tinggi 1403371 1060106   343265
## 3  2022       SMA/SMK/MA 3980426 2670070  1310356
## 4  2022 Perguruan Tinggi 1596304 1211754   384550
## 5  2023       SMA/SMK/MA 3874632 2664972  1209660
## 6  2023 Perguruan Tinggi 1771585 1399198   372387

Interpretasi: Berdasarkan data yang ditampilkan dari tahun 2021 hingga 2023, terlihat bahwa jumlah total lulusan dan jumlah yang termasuk dalam kategori Angkatan Kerja (AK) meningkat untuk kedua jenjang pendidikan, yakni SMA/SMK/MA dan Perguruan Tinggi. Namun, persentase AK dari total lulusan cenderung lebih tinggi pada lulusan Perguruan Tinggi dibandingkan lulusan SMA/SMK/MA setiap tahunnya. Misalnya, pada tahun 2023, 78,98% lulusan Perguruan Tinggi masuk kategori AK, sedangkan untuk lulusan SMA/SMK/MA hanya 68,78%. Hal ini menunjukkan bahwa semakin tinggi jenjang pendidikan terakhir seseorang, semakin besar kemungkinan mereka tergolong sebagai Angkatan Kerja. Selain itu, jumlah Bukan AK cenderung lebih tinggi pada lulusan SMA/SMK/MA, yang bisa mengindikasikan bahwa lulusan jenjang ini memiliki hambatan lebih besar dalam masuk ke pasar kerja dibanding lulusan Perguruan Tinggi.

Tabel Parsial dan Marginal

##           Status
## Pendidikan      AK Bukan AK
##        PT  1060106  2569676
##        SMA  343265  1268526

Interpretasi: Berdasarkan tabel tersebut, terlihat bahwa lulusan Perguruan Tinggi (PT) yang termasuk dalam kategori Angkatan Kerja (AK) berjumlah 1.060.106 orang, jauh lebih tinggi dibanding lulusan SMA yang hanya 343.265 orang. Sebaliknya, jumlah lulusan SMA yang termasuk dalam kategori Bukan AK (1.268.526 orang) jauh lebih besar dibanding lulusan PT yang bukan AK (256.9676 orang). Hal ini menunjukkan bahwa lulusan Perguruan Tinggi memiliki kemungkinan yang lebih besar untuk masuk ke dalam kategori Angkatan Kerja dibanding lulusan SMA. Secara tidak langsung, hal ini juga mengindikasikan bahwa semakin tinggi jenjang pendidikan, semakin tinggi pula keterlibatan seseorang dalam dunia kerja atau aktivitas produktif lainnya.

Tabel Frekuensi Marginal

## $Marginal_Pendidikan
##       PT      SMA 
## 11575776  4888744 
## 
## $Marginal_Tahun
##    2021    2022    2023 
## 5241573 5576730 5646217

Interpretasi: Berdasarkan nilai marginal pada tabel tersebut, dapat dilihat bahwa secara keseluruhan jumlah individu dengan latar belakang pendidikan Perguruan Tinggi (PT) jauh lebih banyak (11.575.776 orang) dibandingkan dengan lulusan SMA/SMK/MA (4.888.744 orang) selama periode waktu yang diamati. Jika ditinjau berdasarkan tahun, jumlah total data juga mengalami peningkatan dari tahun ke tahun: sebanyak 5.241.573 individu pada tahun 2021, meningkat menjadi 5.576.730 pada tahun 2022, dan kembali naik menjadi 5.646.217 pada tahun 2023. Ini menunjukkan adanya tren kenaikan jumlah individu baik dari sisi waktu maupun dari segi jenjang pendidikan, dengan dominasi jumlah dari kelompok lulusan Perguruan Tinggi.

7.2 Distribusi Peluang

Peluang Bersama

##           Status
## Pendidikan         AK   Bukan AK
##        PT  0.06438730 0.15607354
##        SMA 0.02084877 0.07704604

Interpretasi: Tabel tersebut menunjukkan distribusi proporsi gabungan (joint probability) antara tingkat pendidikan dan status pekerjaan (AK dan Bukan AK). Dari hasil ini, terlihat bahwa proporsi tertinggi berasal dari kelompok Perguruan Tinggi yang Bukan AK (0.156), disusul oleh PT yang AK (0.064), SMA yang Bukan AK (0.077), dan terakhir SMA yang AK (0.021). Artinya, meskipun lulusan PT lebih besar kemungkinan menjadi angkatan kerja dibandingkan SMA, jumlah lulusan PT yang belum bekerja atau tidak termasuk AK juga cukup besar secara proporsi. Ini bisa mengindikasikan bahwa meskipun jenjang pendidikan lebih tinggi cenderung meningkatkan peluang masuk angkatan kerja, masih terdapat porsi signifikan lulusan PT yang belum bekerja atau belum termasuk AK.

Tabel Peluang Bersama

##                     Tahun       2021       2022       2023
## Pendidikan Status                                         
## PT         AK             0.06438730 0.07359789 0.08498262
##            Bukan AK       0.15607354 0.16217114 0.16186151
## SMA        AK             0.02084877 0.02335628 0.02261754
##            Bukan AK       0.07704604 0.07958665 0.07347071

Interpretasi: Tabel tersebut menunjukkan peluang gabungan antara variabel Pendidikan, Status (AK atau Bukan AK), dan Tahun (2021–2023). Dari data ini terlihat bahwa proporsi lulusan Perguruan Tinggi yang termasuk AK meningkat tiap tahun, dari 0.064 pada 2021 menjadi 0.085 pada 2023, yang menunjukkan tren positif dalam partisipasi angkatan kerja di kelompok ini. Namun, proporsi Perguruan Tinggi yang Bukan AK masih tetap tinggi, meskipun relatif stabil, sekitar 0.156–0.162. Sementara itu, lulusan SMA/SMK/MA yang masuk AK juga menunjukkan kenaikan kecil pada 2022 (0.023), namun menurun kembali pada 2023 (0.022). Adapun proporsi Bukan AK dari lulusan SMA menurun dari 0.077 di 2021 menjadi 0.073 di 2023, meskipun penurunannya kecil. Secara keseluruhan, data ini mengindikasikan bahwa peningkatan partisipasi AK cenderung lebih tajam pada lulusan PT dibandingkan SMA, namun jumlah Bukan AK di kedua jenjang masih signifikan.

Peluang Marginal

## $P_AK
## [1] 0.6521
## 
## $P_Bukan_AK
## [1] 0.3479

Interpretasi: Nilai probabilitas yang ditampilkan menunjukkan bahwa secara keseluruhan, peluang seseorang termasuk dalam kategori Angkatan Kerja (AK) adalah sebesar 65,21%, sedangkan peluang untuk tidak termasuk dalam Angkatan Kerja (Bukan AK) adalah sebesar 34,79%. Ini berarti mayoritas individu dalam populasi yang diamati tergolong sebagai angkatan kerja. Selisih yang cukup besar antara keduanya juga mencerminkan bahwa partisipasi kerja cukup tinggi dalam data tersebut, namun masih ada sepertiga populasi yang belum masuk dalam kategori angkatan kerja, yang bisa jadi mencakup pelajar, ibu rumah tangga, lansia, atau pengangguran.

Tabel Peluang Marginal

## $Marginal_Pendidikan
##     PT    SMA 
## 0.7031 0.2969 
## 
## $Marginal_Tahun
##   2021   2022   2023 
## 0.3184 0.3387 0.3429

Interpretasi: Berdasarkan hasil marginal probabilitas yang ditampilkan, dapat disimpulkan bahwa mayoritas populasi dalam data berasal dari jenjang Pendidikan Tinggi (PT) dengan proporsi sebesar 70,31%, sementara sisanya 29,69% berasal dari jenjang SMA/SMK/MA. Jika dilihat dari aspek waktu, distribusi tahun relatif seimbang, dengan 34,29% data berasal dari tahun 2023, disusul 33,87% dari tahun 2022, dan 31,84% dari tahun 2021. Ini menunjukkan tren yang stabil dan cukup merata dalam distribusi data lintas waktu, namun secara keseluruhan, pendidikan tinggi mendominasi dalam populasi yang diamati.

Peluang Bersyarat

contoh P(Status = AK l Pendidikan = PT)

## [1] 0.3171328

Interpretasi: Nilai probabilitas bersyarat P(AK | PT) = 0.3171 menunjukkan bahwa di antara seluruh individu yang berasal dari jenjang Pendidikan Tinggi (PT), sekitar 31,71% di antaranya berada dalam status Angkatan Kerja (AK). Ini mengindikasikan bahwa meskipun PT merupakan kelompok pendidikan dengan proporsi terbesar, tidak semua lulusannya langsung masuk ke dalam angkatan kerja, dan masih ada sebagian besar yang berada di luar kategori AK (misalnya belum bekerja, melanjutkan studi, atau tidak aktif bekerja).

P(Tahun = 2023 l AK)

## [1] 0.3713034

Interpretasi: Nilai P(2023 | AK) = 0.3713 berarti bahwa sekitar 37,13% dari seluruh individu yang termasuk dalam Angkatan Kerja (AK) berasal dari tahun 2023. Ini menunjukkan bahwa proporsi angkatan kerja terbanyak berasal dari tahun terakhir dalam data, yang bisa mengindikasikan peningkatan partisipasi angkatan kerja atau jumlah lulusan yang lebih besar pada tahun tersebut dibandingkan tahun-tahun sebelumnya.

Interpretasi: Berdasarkan hasil perhitungan ukuran asosiasi, diperoleh Risk Difference (RD) sebesar 0.102, yang menunjukkan bahwa selisih proporsi individu yang bekerja (AK) antara kelompok pendidikan tinggi dan pendidikan menengah adalah 10,2%. Artinya, individu dengan pendidikan tinggi memiliki peluang lebih tinggi sebesar 10,2% untuk masuk dalam angkatan kerja dibandingkan mereka yang hanya berpendidikan menengah. Relative Risk (RR) sebesar 1.1483 mengindikasikan bahwa peluang individu dari pendidikan tinggi untuk bekerja 1,15 kali lebih besar dibandingkan yang berpendidikan menengah. Sedangkan Odds Ratio (OR) sebesar 1.7055 menunjukkan bahwa perbandingan odds bekerja antara kelompok pendidikan tinggi dan menengah adalah 1,71 banding 1, sehingga dapat disimpulkan bahwa pendidikan tinggi secara signifikan meningkatkan peluang individu untuk menjadi bagian dari angkatan kerja.

Struktur Tabel Kontingensi 3D (AK × Pendidikan × Tahun)

# Asumsikan data_ak sudah dibaca dan diolah seperti ini:
# (JANGAN JALANKAN KEMBALI JIKA SUDAH ADA)
# data_ak <- read.csv("D:/SEMESTER 4/ADK/Data Pendidikan Terakhir.csv", skip = 2)  
# data_ak <- data_ak[1:4, ]
# colnames(data_ak)[c(1, 4, 7)] <- c("Pendidikan", "Total_2023", "Persen_AK_2023")
# data_ak$Total_2023 <- as.numeric(as.character(data_ak$Total_2023))
# data_ak$Persen_AK_2023 <- as.numeric(as.character(data_ak$Persen_AK_2023))
# data_ak$AK <- round(data_ak$Total_2023 * data_ak$Persen_AK_2023 / 100)
# data_ak$Bukan_AK <- data_ak$Total_2023 - data_ak$AK

# Buat array 3D: [Status AK, Pendidikan, Tahun]
data_array <- array(
  data = c(data_ak$AK, data_ak$Bukan_AK),
  dim = c(2, 4, 1),
  dimnames = list(
    Status = c("AK", "Bukan AK"),
    Pendidikan = data_ak$Pendidikan,
    Tahun = "2023"
  )
)

data_array
## , , Tahun = 2023
## 
##           Pendidikan
## Status     >= SD/MI SMP/MTs SMA/SMK/MA Perguruan Tinggi
##   AK         638740 2664972     505955          1209660
##   Bukan AK   724330 1399198     807671           372387
# Tabel marginal (menjumlahkan sepanjang dimensi tahun)
tabel_marginal <- apply(data_array, c(1, 2), sum)

# Tampilkan tabel dalam bentuk data frame
tabel_marginal_df <- as.data.frame(t(tabel_marginal))
tabel_marginal_df
# PT = baris ke-4, SMA = baris ke-3
n11 <- data_ak$AK[4]        # PT - AK
n12 <- data_ak$Bukan_AK[4]  # PT - Bukan AK
n21 <- data_ak$AK[3]        # SMA - AK
n22 <- data_ak$Bukan_AK[3]  # SMA - Bukan AK

# Proporsi
p1 <- n11 / (n11 + n12)
p2 <- n21 / (n21 + n22)

# Ukuran asosiasi
RD <- p1 - p2
RR <- p1 / p2
OR <- (n11 / n12) / (n21 / n22)

# Buat list hasil
hasil <- list(
  Risk_Difference = round(RD, 4),
  Relative_Risk = round(RR, 4),
  Odds_Ratio = round(OR, 4)
)

# Tampilkan
hasil
## $Risk_Difference
## [1] 0.102
## 
## $Relative_Risk
## [1] 1.1483
## 
## $Odds_Ratio
## [1] 1.7055

Interpretasi: Berdasarkan data pendidikan terakhir tahun 2023 dan hasil pengukuran asosiasi, terlihat bahwa individu dengan pendidikan terakhir perguruan tinggi memiliki peluang lebih tinggi untuk masuk kategori angkatan kerja (AK) dibandingkan dengan mereka yang hanya lulusan SMA/SMK/MA. Nilai Risk Difference (RD) sebesar 0.102 menunjukkan bahwa selisih proporsi AK antara kelompok pendidikan tinggi dan SMA/SMK adalah 10,2%. Relative Risk (RR) sebesar 1.1483 menunjukkan bahwa lulusan perguruan tinggi memiliki risiko 14,83% lebih besar untuk masuk ke dalam AK dibandingkan dengan lulusan SMA/SMK. Sementara itu, Odds Ratio (OR) sebesar 1.7055 memperlihatkan bahwa odds (kemungkinan relatif) menjadi AK bagi lulusan perguruan tinggi sekitar 1,7 kali lipat lebih besar dibandingkan lulusan SMA/SMK. Hasil ini menunjukkan adanya asosiasi positif antara tingkat pendidikan dan partisipasi dalam angkatan kerja, di mana semakin tinggi pendidikan, semakin besar kemungkinan seseorang tergolong ke dalam AK.

data3_cmh <- array(
  data = c(
    # 2021
    1060106, 2569676,
    343265,  1268526,
    # 2022
    1211754, 2670070,
    384550,  1310356,
    # 2023
    1399198, 2664972,
    372387,  1209660
  ),
  dim = c(2, 2, 3),
  dimnames = list(
    Pendidikan = c("PT", "SMA"),
    Status = c("AK", "Bukan AK"),
    Tahun = c("2021", "2022", "2023")
  )
)

7. 3 Ukuran Asosiasi (RD, RR, OR)

for (tahun in dimnames(data3_cmh)$Tahun) {
  a <- data3_cmh["PT", "AK", tahun]
  b <- data3_cmh["PT", "Bukan AK", tahun]
  c <- data3_cmh["SMA", "AK", tahun]
  d <- data3_cmh["SMA", "Bukan AK", tahun]
  
  p1 <- a / (a + b)
  p2 <- c / (c + d)
  
  RD <- p1 - p2
  RR <- p1 / p2
  OR <- (a * d) / (b * c)
  
  cat(paste0("\nTahun ", tahun, "\n"))
  cat("Risk Difference (RD): ", round(RD, 4), "\n")
  cat("Relative Risk (RR): ", round(RR, 4), "\n")
  cat("Odds Ratio (OR): ", round(OR, 4), "\n")
}
## 
## Tahun 2021
## Risk Difference (RD):  0.0859 
## Relative Risk (RR):  1.1283 
## Odds Ratio (OR):  1.5245 
## 
## Tahun 2022
## Risk Difference (RD):  0.0883 
## Relative Risk (RR):  1.1316 
## Odds Ratio (OR):  1.5464 
## 
## Tahun 2023
## Risk Difference (RD):  0.102 
## Relative Risk (RR):  1.1483 
## Odds Ratio (OR):  1.7055

Interpretasi: Pada tahun 2023, Risk Difference (RD) menunjukkan bahwa proporsi angkatan kerja lulusan PT lebih tinggi sekitar 10.2% dibandingkan lulusan SMA. Relative Risk (RR) sebesar 1.1483 berarti lulusan PT memiliki peluang 14.8% lebih besar untuk masuk angkatan kerja dibanding lulusan SMA. Odds Ratio (OR) sebesar 1.7055 menunjukkan odds lulusan PT masuk angkatan kerja 1.7 kali lebih besar dari lulusan SMA.

7.4 CMH

cmh_result <- mantelhaen.test(data3_cmh, correct = FALSE)
cmh_result
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  data3_cmh
## Mantel-Haenszel X-squared = 138847, df = 1, p-value < 2.2e-16
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.588484 1.596316
## sample estimates:
## common odds ratio 
##          1.592395

Interpretasi: Karena p-value dari uji Cochran-Mantel-Haenszel (CMH) < 0.05, maka kita menolak H₀, yang berarti ada hubungan signifikan antara pendidikan dan status angkatan kerja setelah mengontrol tahun. Artinya, meskipun waktu (tahun) diperhitungkan, tingkat pendidikan tetap berpengaruh signifikan terhadap status AK seseorang.

8 Generalized Linear Model (GLM)

Generalized Linear Model (GLM) merupakan perluasan dari model regresi linear klasik. GLM memungkinkan digunakan untuk memodelkan data dimana variabel respons tidak berdistribusi normal dan atau hubungan antara variabel prediktor dengan rata-rata dari variabel respons tidak linear. GLM terdiri dari tiga komponen utama:

  1. Distribusi dari exponential family untuk variabel respons.
  2. Fungsi link yang menghubungkan ekspektasi dari variabel respons ke kombinasi linear dari variabel prediktor.
  3. Fungsi linear prediktor : 𝜂 =𝑋𝛽

8.1 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 yang Termasuk Exponential Family:

  • Distribusi Normal
  • Distribusi Binomial
  • Distribusi Poisson
  • Distribusi Gamma

Contoh Pembuktian: Distribusi Binomial

Fungsi probabilitas distribusi binomial adalah:

\[ P(Y = y) = \binom{n}{y} \pi^y (1 - \pi)^{n - y} \]

Tuliskan ulang dalam bentuk exponential family:

\[ P(Y = y) = \exp \left\{ \log \binom{n}{y} + y \log \left( \frac{\pi}{1 - \pi} \right) + n \log (1 - \pi) \right\} \]

Dengan substitusi:

  • \(\theta = \log \left( \frac{\pi}{1 - \pi} \right)\)
  • \(b(\theta) = -n \log(1 - \pi)\)
  • \(\phi = 1\)

Kesimpulan : Dengan bentuk tersebut, maka distribusi binomial termasuk dalam keluarga exponential family.

8.2 Model Regresi Logistik

Regresi Logistik memiliki kemiripan dengan regresi linear, yaitu menggabungkan nilai-nilai input secara linear dengan bobot (koefisien) untuk menghasilkan sebuah prediksi. Bedanya, regresi logistik membatasi hasil prediksinya agar berada dalam bentuk nilai biner 0 atau 1 dengan bantuan fungsi sigmoid sebagai fungsi aktivasi. Fungsi ini mengubah output menjadi nilai antara 0 sampai 1 dan membentuk kurva seperti huruf S.

Model ini digunakan untuk menganalisis hubungan antara satu atau lebih variabel independen dan mengklasifikasikan data ke dalam kelas-kelas yang bersifat diskrit. Artinya, regresi logistik cocok digunakan untuk masalah klasifikasi, terutama yang hanya memiliki dua kemungkinan kategori (dikenal dengan klasifikasi biner).

Contohnya, angka 0 bisa merepresentasikan gagal, sedangkan angka 1 merepresentasikan sukses.

Beberapa contoh penerapan regresi logistik dalam klasifikasi biner antara lain:

-Memprediksi kelulusan siswa: Model regresi logistik dapat digunakan untuk memperkirakan apakah seorang siswa akan lulus atau tidak berdasarkan variabel-variabel seperti kehadiran, nilai tugas, partisipasi kelas, dan nilai ujian.

-Deteksi penipuan transaksi keuangan: Dalam dunia perbankan dan keuangan, regresi logistik dapat membantu mendeteksi apakah sebuah transaksi tergolong normal atau mencurigakan (fraud) dengan menganalisis pola transaksi seperti frekuensi, jumlah uang, lokasi transaksi, dan waktu.

-Memprediksi hasil diagnosis penyakit: Dalam dunia medis, model ini bisa digunakan untuk menentukan apakah pasien menderita suatu penyakit (misalnya diabetes atau kanker) berdasarkan variabel seperti usia, tekanan darah, kadar gula, dan riwayat keluarga.

Keunggulan Utama Regresi Logistik yaitu mudah diterapkan dalam Machine Learning. Regresi logistik sangat cocok digunakan dalam Machine Learning karena proses pelatihan (training) dan pengujian (testing) nya cukup sederhana. Model ini belajar dari pola-pola dalam data input dan menghubungkannya dengan output (label). Karena tidak membutuhkan komputasi tinggi, regresi logistik tergolong mudah untuk diimplementasikan, dipahami, dan dilatih dibanding metode lain dalam Machine Learning.

Cocok untuk data yang bisa dipisahkan secara linear Jika dua kelas data bisa dipisahkan dengan garis lurus (linear), maka regresi logistik akan sangat efektif dalam mengklasifikasikan data tersebut ke dalam dua kelompok berbeda. Dalam konteks ini, variabel target (y) hanya memiliki dua nilai, sehingga model bisa menentukan batas yang jelas di antara keduanya.

Memberikan wawasan yang bermakna Regresi logistik juga bisa menunjukkan seberapa besar pengaruh suatu variabel terhadap hasil akhir, melalui koefisiennya. Selain itu, model ini juga menunjukkan apakah hubungan antar variabel bersifat positif atau negatif.

Persamaan dan Asumsi dalam Regresi Logistik Model ini menggunakan fungsi sigmoid, yaitu fungsi matematika berbentuk kurva S yang mengubah nilai output ke dalam rentang antara 0 dan 1. Jika hasil dari fungsi sigmoid (yang merepresentasikan probabilitas) lebih tinggi dari ambang batas tertentu, maka model akan mengklasifikasikan data sebagai bagian dari suatu kelas. Sebaliknya, jika di bawah ambang tersebut, maka data dianggap tidak termasuk dalam kelas itu.

Selain itu, jika keluaran dari fungsi sigmoid (yaitu probabilitas yang diperkirakan) lebih besar dari ambang batas yang telah ditentukan dalam grafik, maka model akan memprediksi bahwa suatu observasi termasuk dalam kelas tertentu. Sebaliknya, jika nilai probabilitas tersebut lebih kecil dari ambang batas, model akan memprediksi bahwa observasi tersebut tidak termasuk ke dalam kelas tersebut

Sebagai contoh:

-Jika hasil fungsi sigmoid lebih dari 0,5, maka output dianggap sebagai 1 (kelas positif). -Jika hasilnya kurang dari 0,5, maka output diklasifikasikan sebagai 0 (kelas negatif). -Jika grafik menuju ke arah negatif secara ekstrem, maka nilai prediksi y akan menjadi 0,dan sebaliknya.

Fungsi sigmoid digunakan dalam regresi logistik untuk mengubah nilai prediksi menjadi probabilitas. Rumus dari fungsi sigmoid yaitu: \[ f(x) = \frac{1}{1 + e^{-x}} \] Persiapan Data

# Data Pendidikan dan Status AK
pendidikan <- c("SD/MI", "SMP/MTs", "SMA/SMK/MA", "Perguruan Tinggi")
total <- c(1144695, 1532001, 3874632, 1771585)
persen_ak <- c(55.76, 46.86, 66.95, 75.54)

# Hitung AK dan Bukan AK
data_ak <- data.frame(
  Pendidikan = pendidikan,
  Total = total,
  AK = round(total * persen_ak / 100)
)
data_ak$Bukan_AK <- data_ak$Total - data_ak$AK

# Konversi jadi data biner (per observasi)
data_long <- data.frame(
  Pendidikan = rep(data_ak$Pendidikan, times = data_ak$Total),
  Status = unlist(mapply(function(ak, bukan) c(rep(1, ak), rep(0, bukan)), 
                         data_ak$AK, data_ak$Bukan_AK))
)

8.2.1 Regresi Logistik

Kurva sigmoid dalam regresi logistik menunjukkan hubungan non-linear antaravariabel prediktor dan probabilitas output. Pendekatan ini efektif untuk klasifikasi biner seperti deteksi penyakit, kelulusan siswa, dan prediksi ya/tidak.

Fungsi Link (Logit): Fungsi link logit dapat dinyatakan dalam bentuk berikut:

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

Model Regresi Logistik:

\[ \log \left( \frac{\mu}{1 - \mu} \right) = 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_{i=1}^{n} \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right] \] Dengan:

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

# Model GLM binomial (logistik)
model_logistik <- glm(Status ~ Pendidikan, data = data_long, family = binomial())
summary(model_logistik)
## 
## Call:
## glm(formula = Status ~ Pendidikan, family = binomial(), data = data_long)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           1.127622   0.001748   645.2   <2e-16 ***
## PendidikanSD/MI      -0.896195   0.002568  -348.9   <2e-16 ***
## PendidikanSMA/SMK/MA -0.421698   0.002055  -205.2   <2e-16 ***
## PendidikanSMP/MTs    -1.253387   0.002382  -526.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 10919845  on 8322912  degrees of freedom
## Residual deviance: 10577728  on 8322909  degrees of freedom
## AIC: 10577736
## 
## Number of Fisher Scoring iterations: 4
# Transformasi ke Odds Ratio
exp(coef(model_logistik))
##          (Intercept)      PendidikanSD/MI PendidikanSMA/SMK/MA 
##            3.0883045            0.4081198            0.6559322 
##    PendidikanSMP/MTs 
##            0.2855361

Interpretasi Koefisien

• Intercept ((Intercept)): nilai log-odds saat x=0

• Koefisien x: perubahan log-odds untuk setiap satu satuan peningkatan pada x

• Nilai p: menunjukkan signifikansi statistik dari prediktor

Visualisasi Hasil Model Logistik

menggunakan ggplot dengan syntax: ggplot(data_long, aes(x = Pendidikan, y = Status)) + geom_jitter(width = 0.2, height = 0.05, alpha = 0.4) + stat_summary(fun = mean, geom = “point”, aes(y = pred), color = “blue”, size = 3) + labs(title = “Prediksi Probabilitas AK per Jenjang Pendidikan”, y = “Status AK (1 = AK, 0 = Bukan AK)”) + theme_minimal()

8.3 Regresi Poisson (Model Jumlah AK)

Regresi Poisson digunakan ketika variabel respons adalah data cacah (count data),yaitu bilangan bulat non negatif. Model ini merupakan bagian dari Generalized Linear Model(GLM) dengan asumsi bahwa distribusi variabel respons adalah distribusi Poisson.

Distribusi Poisson memiliki fungsi probabilitas:

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

Distribusi Poisson dapat dituliskan dalam bentuk exponential family sebagai berikut:

\[ f(y; \theta) = \exp \left\{ y \log(\lambda) - \lambda - \log(y!) \right\} \]

dengan:

  • \[ \theta = \log(\lambda) \]
  • \[ b(\theta) = e^{\theta} = \lambda \]
  • \[ \phi = 1 \]
  • \[ c(y, \phi) = -\log(y!) \]

Maka distribusi Poisson termasuk dalam exponential family.

Fungsi Link

Fungsi link kanonik untuk distribusi Poisson adalah fungsi logaritma:

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

Sehingga modelnya menjadi:

\[ \log(\mu_i) = x_i^T \beta \]

dan fungsi inverse link:

\[ \mu_i = \exp(x_i^T \beta) \]

Estimasi Parameter

Estimasi parameter \(( \beta )\) dilakukan dengan metode Maximum Likelihood Estimation (MLE).
Log-likelihood fungsi untuk regresi Poisson:

\[ l(\beta) = \sum_{i=1}^n \left[ y_i x_i^T \beta - \exp(x_i^T \beta) - \log(y_i!) \right] \]

Nilai \(( \beta )\) dapat diperoleh melalui metode numerik seperti iterasi Newton-Raphson.

model_poisson <- glm(AK ~ Pendidikan, data = data_ak, family = poisson())
summary(model_poisson)
## 
## Call:
## glm(formula = AK ~ Pendidikan, family = poisson(), data = data_ak)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          14.1068771  0.0008644 16319.3   <2e-16 ***
## PendidikanSD/MI      -0.7403416  0.0015212  -486.7   <2e-16 ***
## PendidikanSMA/SMK/MA  0.6618600  0.0010643   621.9   <2e-16 ***
## PendidikanSMP/MTs    -0.6227971  0.0014629  -425.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1.7228e+06  on 3  degrees of freedom
## Residual deviance: 4.0824e-10  on 0  degrees of freedom
## AIC: 71.078
## 
## Number of Fisher Scoring iterations: 2
# Interpretasi Koefisien Poisson
exp(coef(model_poisson))
##          (Intercept)      PendidikanSD/MI PendidikanSMA/SMK/MA 
##         1.338255e+06         4.769510e-01         1.938394e+00 
##    PendidikanSMP/MTs 
##         5.364419e-01

8.4 Diagnostik Model Poisson

Salah satu asumsi penting dari model Poisson adalah bahwa mean dan varians dari variabel respons adalah sama (𝐸[𝑌]=𝑉𝑎𝑟[𝑌]). Jika varians lebih besar dari mean, maka terjadi overdispersion. Untuk mendeteksi overdispersion:

dispersion <- sum(residuals(model_poisson, type = "pearson")^2) / model_poisson$df.residual
dispersion
## [1] Inf

Jika nilai dispersion > 1, maka overdispersion mungkin terjadi dan model alternatif seperti Negative Binomial Regression dapat digunakan. visualisasi: pakai syntax ini ggplot(data_ak, aes(x = Pendidikan, y = AK)) + geom_col(fill = “gray”) + geom_point(aes(y = Prediksi), color = “red”, size = 3) + labs(title = “Prediksi Jumlah AK per Pendidikan (Model Poisson)”, y = “Jumlah AK”) + theme_minimal()

8.5 Kesimpulan

Model GLM (logistik dan Poisson) dapat digunakan untuk menganalisis hubungan antara jenjang pendidikan dan status angkatan kerja. Regresi logistik cocok untuk memodelkan probabilitas individu masuk AK, sedangkan regresi Poisson cocok untuk jumlah AK pada masing-masing kategori pendidikan.

9 Inferensi GLM

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

Ekspektasi dan Varians dalam GLM

  1. Ekspektasi Estimator

Ekspektasi menunjukkan apakah suatu estimator tak bias, yaitu:

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

Dalam GLM, MLE dari \(\hat{\beta}\) bersifat asymptotically unbiased.

  1. Varians Estimator

Varians menunjukkan presisi dari estimasi parameter:

\[ \text{Var}(\hat{\beta}) \approx (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \]

dimana \(\mathbf{W}\) adalah matriks bobot yang tergantung pada distribusi dan fungsi link.

Distribusi Asimptotik Estimator

Dengan ukuran sampel besar:

\[ \hat{\beta} \sim \mathcal{N}(\beta, \text{Var}(\hat{\beta})) \]

Distribusi ini adalah dasar dari:

Varians dalam GLM Tidak Konstan

Tidak seperti regresi linear (OLS) yang mengasumsikan homoskedastisitas:

\[ \text{Var}(Y_i) = \sigma^2 \]

Dalam GLM:

\[ \text{Var}(Y_i) = \phi V(\mu_i) \]

Contoh:

Perhitungan menggunakan R

Kesimpulan:

9.1 Mencari Ekspektasi dan Varians dalam GLM

Ekspektasi

Jika diturunkan berdasarkan fungsi momen:

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

Untuk keluarga eksponensial:

\[ \log f(y; \theta) = a(y) + b(\theta) y + c(\theta) \]

atau:

\[ \log f(y; \theta) = y \theta - b(\theta) + c(y) \]

Maka ekspektasi turunan pertama:

\[ U(\theta) = \frac{\partial \ell}{\partial \theta} = y - b'(\theta) \] Dan ekspektasi turunan pertama:

\[ \mathbb{E}[U(\theta)] = \mathbb{E}[y - b'(\theta)] = \mu - b'(\theta) = 0 \]

Maka:

\[ \mu = b'(\theta) \]

Varians
Turunan kedua:

\[ \frac{\partial^2 \ell}{\partial \theta^2} = -b''(\theta) \]

Sehingga:

\[ \text{Var}(Y) = b''(\theta) = \phi V(\mu) \]

9.2 Metode Penaksiran Parameter

Maximum Likelihood Estimation (MLE)

  • Prinsip dasar: memaksimumkan fungsi likelihood/log-likelihood.
  • Langkah:
    • Turunan pertama = 0
    • Turunan kedua < 0

Namun, karena bentuk GLM tidak eksplisit, digunakan metode numerik.

Metode Optimisasi: Newton-Raphson

  • Menggunakan score vector (gradien)
  • Menggunakan Hessian matrix

Iterasi:

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

  • Modifikasi Newton-Raphson, mengganti Hessian dengan matriks informasi Fisher.

IRLS (Iteratively Reweighted Least Square) - Modifikasi dari Fisher scoring, hasil estimasi mirip dengan Least Square.

Implementasi Newton-Raphson Statistik score ke-𝑗 Statistik score ke-𝑗:

\[ U_j(\beta) = \frac{\partial \log L(\beta)}{\partial \beta_j} \] Turunan kedua:

\[ H_{jk}(\beta) = \frac{\partial^2 \log L(\beta)}{\partial \beta_j \partial \beta_k} \]

Taylor expansion:

\[ U(\beta^*) \approx U(\beta) + H(\beta)(\beta^* - \beta) \]

Estimasi parameter:

\[ \hat{\beta} \approx \beta^{(t)} - H^{-1}(\beta^{(t)}) U(\beta^{(t)}) \]

9.3 Diagnostik Model GLM

Diagnostik digunakan untuk mengevaluasi apakah model sudah tepat.

  • Uji formal
  • Grafik antara nilai prediksi vs nilai aktual

Statistik Devians

  • Mengukur apakah ada model lain yang lebih baik.
  • Nilai devians besar → model tidak cocok.
  • Devians adalah:

\[ D = 2 \sum \left[ y_i \log \left( \frac{y_i}{\hat{\mu}_i} \right) - (y_i - \hat{\mu}_i) \right] \]

  • Devians membandingkan model terhadap saturated model.
  • Devians kecil → model lebih cocok.

Statistik Chi-Kuadrat Pearson

  • Menguji apakah model lebih baik daripada tidak ada model sama sekali.
  • Statistik:

\[ X^2 = \sum \left( \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i} \right) \] • Jika signifikan → model lebih baik daripada tanpa model.

Catatan

• Untuk data yang dikelompokkan, statistik devians dan chi-kuadrat Pearson mengikuti distribusi ChiSquare.

• Untuk data tidak dikelompokkan, tidak mengikuti distribusi Chi-Square.

• Devians diminimalkan oleh MLE → cocok digunakan untuk evaluasi model.

Analisis Residual

• Residual adalah selisih antara observasi dengan prediksi.

• Dapat digunakan untuk memeriksa penyimpangan sistematis.

• Dapat diplot untuk menilai asumsi model

9.4 Detail Metode Estimasi dan Inferensi Regresi Logistik

Regresi logistik digunakan untuk memodelkan probabilitas dari variabel respons biner (0/1) berdasarkan satu atau lebih variabel prediktor. Estimasi parameter dilakukan menggunakan Maximum Likelihood Estimation (MLE) karena model tidak linear dalam parameternya.

Fungsi model logistik:

\[ \pi(x) = \frac{\exp(\beta_0 + \beta_1 x)}{1 + \exp(\beta_0 + \beta_1 x)} \]

Log-likelihood untuk \(( n )\) observasi:

\[ \ell(\beta) = \sum_{i=1}^{n} \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right] \] Estimasi dengan Newton-Raphson

Metode Newton-Raphson digunakan untuk mencari nilai parameter \(( \beta )\) yang memaksimalkan fungsi log-likelihood pada model regresi logistik.

Fungsi Log-Likelihood

Model regresi logistik untuk probabilitas:

\[ \pi_i = \frac{1}{1 + \exp(-x_i^\top \beta)} \]

Log-likelihood untuk \(( n )\) observasi:

\[ \ell(\beta) = \sum_{i=1}^{n} \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right] \]

Langkah-langkah Newton-Raphson Turunan dalam Estimasi Logistik

  1. Turunan Pertama (Score Function)

\[ U(\beta) = \frac{\partial \ell(\beta)}{\partial \beta} = X^\top (y - \pi) \]

  1. Turunan Kedua (Hessian Matrix)

\[ H(\beta) = -X^\top W X, \quad \text{dengan } W = \text{diag}(\pi_i (1 - \pi_i)) \]

  1. Iterasi Newton-Raphson

\[ \beta^{(t+1)} = \beta^{(t)} + (X^\top W^{(t)} X)^{-1} X^\top (y - \pi^{(t)}) \]

Perhitungan menggunakan R

# Data
data <- data.frame(
  Pendidikan = c("SD", "SMP", "SMA", "PT"),
  AK = c(638740, 724330, 2664972, 1399198),
  Total = c(1144695, 1532001, 3874632, 1771585)
)

# Tambahkan log Total sebagai offset (optional jika dibutuhkan)
data$log_total <- log(data$Total)

# Ubah Pendidikan jadi faktor
data$Pendidikan <- factor(data$Pendidikan)

# Regresi Poisson: respons = AK, prediktor = Pendidikan
model_poisson <- glm(AK ~ Pendidikan, family = poisson(), data = data)

# Ringkasan model
summary(model_poisson)
## 
## Call:
## glm(formula = AK ~ Pendidikan, family = poisson(), data = data)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   14.1514098  0.0008454 16739.4   <2e-16 ***
## PendidikanSD  -0.7841570  0.0015101  -519.3   <2e-16 ***
## PendidikanSMA  0.6442943  0.0010440   617.1   <2e-16 ***
## PendidikanSMP -0.6584074  0.0014475  -454.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  1.8124e+06  on 3  degrees of freedom
## Residual deviance: -5.0906e-10  on 0  degrees of freedom
## AIC: 71.159
## 
## Number of Fisher Scoring iterations: 2

Newton-Raphson Iterasi Manual

# Data
pendidikan <- factor(c("SD", "SMP", "SMA", "PT"), 
                     levels = c("SD", "SMP", "SMA", "PT"), ordered = TRUE)
X1 <- as.numeric(pendidikan)  # representasi ordinal (1,2,3,4)
AK <- c(638740, 724330, 2664972, 1399198)
Bukan_AK <- c(505955, 807671, 1209660, 372387)

# Bentuk data biner
x <- rep(X1, times = AK + Bukan_AK)
y <- c()
for(i in 1:4){
  y <- c(y, rep(1, AK[i]), rep(0, Bukan_AK[i]))
}

# Buat matriks desain
X <- cbind(1, x)  # intercept dan prediktor ordinal pendidikan

# Newton-Raphson Iterasi Manual
beta <- c(0, 0)  # inisialisasi
tol <- 1e-6
max_iter <- 100

Selanjutnya pakai syntax:

for (i in 1:max_iter) { eta <- X %% beta pi_hat <- 1 / (1 + exp(-eta)) W <- diag(as.numeric(pi_hat (1 - pi_hat))) z <- eta + solve(W) %% (y - pi_hat) beta_new <- solve(t(X) %% W %% X) %% t(X) %% W %% z if (sum(abs(beta_new - beta)) < tol) break beta <- beta_new }

beta # hasil estimasi akhir • Estimasi parameter pada model regresi logistik dilakukan dengan MLE. • Newton-Raphson adalah metode numerik yang digunakan untuk memaksimalkan log-likelihood.

• Iterasi didasarkan pada turunan pertama (score) dan kedua (Hessian).

• Prosedur identik dengan IRLS (Iteratively Reweighted Least Squares) dalam implementasi GLM.

Inferensi Parameter

  1. Uji Wald

Tujuan Uji Wald
Untuk menguji signifikansi parameter \(( \beta_j )\) dalam model regresi logistik:

  • \(( H_0: \beta_j = 0 )\)
  • \(( H_1: \beta_j \neq 0 )\)

Teori Wald Test
Dari teori estimasi MLE, estimator \(( \hat{\beta}_j )\) mendekati distribusi normal:

\[ \hat{\beta}_j \sim N(\beta_j, \text{Var}(\hat{\beta}_j)) \]

Jika \(( H_0 )\) benar (yaitu \(( \beta_j = 0 )\)), maka:

\[ Z = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \sim N(0,1) \]

Dengan Statistik Wald:

\[ W = Z^2 = \left( \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \right)^2 \sim \chi^2_1 \]

Uji Wald Langkah demi Langkah

# Data dari Andin
pendidikan <- c("SD/MI", "SMP/MTs", "SMA/SMK/MA", "Perguruan Tinggi")
Total_2023 <- c(1144695, 1532001, 3874632, 1771585)
Persen_AK_2023 <- c(55.76, 46.86, 66.95, 75.54)

# Hitung AK dan Bukan AK
AK <- round(Total_2023 * Persen_AK_2023 / 100)
Bukan_AK <- Total_2023 - AK

# Buat data frame biner
x <- factor(rep(pendidikan, times = AK + Bukan_AK), 
            levels = pendidikan)
y <- unlist(mapply(function(ak, bukan) c(rep(1, ak), rep(0, bukan)), 
                  AK, Bukan_AK))

data <- data.frame(x = x, y = y)
data$x <- relevel(data$x, ref = "SMA/SMK/MA")  # Ref kategori

# Model regresi logistik
model <- glm(y ~ x, data = data, family = binomial)

summary(model)
## 
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.705924   0.001080   653.6   <2e-16 ***
## xSD/MI            -0.474497   0.002170  -218.7   <2e-16 ***
## xSMP/MTs          -0.831689   0.001946  -427.3   <2e-16 ***
## xPerguruan Tinggi  0.421698   0.002055   205.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 10919845  on 8322912  degrees of freedom
## Residual deviance: 10577728  on 8322909  degrees of freedom
## AIC: 10577736
## 
## Number of Fisher Scoring iterations: 4

Langkah 1: Ambil nilai koefisien dan standard error (SE)

beta_hat <- coef(model)["xPerguruan Tinggi"]
se_beta <- summary(model)$coefficients["xPerguruan Tinggi", "Std. Error"]

Langkah 2: Hitung Statistik Z

Z <- beta_hat / se_beta
Z
## xPerguruan Tinggi 
##          205.2466

Langkah 3: Hitung Statistik Wald

Wald_stat <- Z^2
Wald_stat
## xPerguruan Tinggi 
##          42126.16

Langkah 4: Hitung p-value

p_value <- 1 - pchisq(Wald_stat, df = 1)
p_value
## xPerguruan Tinggi 
##                 0

Interpretasi: Interpretasi

• Jika p-value < 0.05, maka koefisien signifikan → variabel prediktor berpengaruh.

• Jika p-value > 0.05, maka tidak ada cukup bukti untuk menolak - \(( H_0 )\)

Kesimpulan

Uji Wald didasarkan pada rasio antara estimasi parameter dan standar error-nya. Dengan menaikkan nilai Z menjadi kuadrat (Z²), kita memperoleh distribusi Chi-Square untuk pengujian hipotesis parameter individual dalam model regresi logistik

  1. Uji Likelihood Ratio (Chi-Square)

Bandingkan model penuh dengan model tanpa prediktor

# Data
data_ak <- data.frame(
  Pendidikan = c(">= SD/MI", "SMP/MTs", "SMA/SMK/MA", "Perguruan Tinggi"),
  Total_2023 = c(1144695, 1532001, 3874632, 1771585),
  Persen_AK_2023 = c(55.76, 46.86, 66.95, 75.54)
)

# Hitung AK dan Bukan AK
data_ak$AK <- round(data_ak$Total_2023 * data_ak$Persen_AK_2023 / 100)
data_ak$Bukan_AK <- data_ak$Total_2023 - data_ak$AK

# Ubah ke format long untuk regresi logistik
data_long <- data_ak %>%
  select(Pendidikan, AK, Bukan_AK) %>%
  tidyr::pivot_longer(cols = c("AK", "Bukan_AK"), names_to = "Status", values_to = "Jumlah")

data_long$Status <- factor(data_long$Status, levels = c("Bukan_AK", "AK"))
data_long$Pendidikan <- factor(data_long$Pendidikan)

# Variabel respons biner
data_long$Y <- ifelse(data_long$Status == "AK", 1, 0)

# Model penuh
model <- glm(Y ~ Pendidikan, data = data_long, weights = Jumlah, family = binomial)

# Model null (tanpa prediktor)
model_null <- glm(Y ~ 1, data = data_long, weights = Jumlah, family = binomial)

# Likelihood ratio test
anova(model_null, model, test = "Chisq")

Evaluasi Kebaikan Model

  1. Akaike Information Criterion (AIC)

Semakin kecil AIC, semakin baik model.

AIC(model)
## [1] 10577736
  1. Bayesian Information Criterion (BIC)

Alternatif terhadap AIC, menghukum kompleksitas model.

BIC(model)
## [1] 10577736

• Estimasi regresi logistik dilakukan dengan MLE melalui iterasi Newton-Raphson.

• Inferensi parameter dapat dilakukan dengan uji Wald dan likelihood ratio (uji Chi-Square).

• AIC dan BIC digunakan untuk mengevaluasi kompleksitas dan kecocokan model.

9.5 Detail Metode Estimasi dan Inferensi Regresi Poisson

Model regresi Poisson digunakan untuk memodelkan data count (jumlah kejadian) dimana variabel respons mengikuti distribusi Poisson. Estimasi dilakukan dengan Maximum Likelihood Estimation (MLE), dan inferensi dilakukan dengan uji Wald dan Likelihood Ratio Test.

Model Regresi Poisson

Distribusi Poisson:

\[ P(Y_i = y_i) = \frac{e^{-\lambda_i} \lambda_i^{y_i}}{y_i!} \]

Model regresi Poisson:

\[ \log(\lambda_i) = x_i^T \beta \] Estimasi Parameter (MLE)

Log-likelihood fungsi:

\[ \ell(\beta) = \sum_{i=1}^{n} \left[ y_i \log(\lambda_i) - \lambda_i - \log(y_i!) \right] \]

Dengan:

\[ \lambda_i = \exp(x_i^T \beta) \]

Estimasi dilakukan dengan metode iterasi (IRLS)

Estimasi parameter model regresi Poisson menggunakan pendekatan
Maximum Likelihood Estimation (MLE) dengan metode
Iteratively Reweighted Least Squares (IRLS) secara manual, tanpa menggunakan glm().

Tahap 1: Definisikan Model Regresi Poisson

\[ \log(\lambda_i) = x_i^T \beta \quad \Rightarrow \quad \lambda_i = \exp(x_i^T \beta) \]

Tahap 2: Mencari Log-likelihood yang dimaksimumkan

\[ \ell(\beta) = \sum_{i=1}^{n} \left[ y_i \log(\lambda_i) - \lambda_i - \log(y_i!) \right] \] Tahap 3: Formulasi iteratif:

\[ \beta^{(t+1)} = \left( \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{z}^{(t)} \]

Dengan:

  • \[\mathbf{W} = \text{diag}(\lambda_i)\]
  • \[\mathbf{z} = \eta + \frac{y - \lambda}{\lambda}\]

dan

\[ \eta_i = \log(\lambda_i) = x_i^T \beta \] Simulasi Data

Perhitungan menggunakan R

# Data AK berdasarkan pendidikan
Pendidikan <- c("SD", "SMP", "SMA", "PT")
Total_2023 <- c(1144695, 1532001, 3874632, 1771585)
Persen_AK_2023 <- c(55.76, 46.86, 66.95, 75.54)

AK <- round(Total_2023 * Persen_AK_2023 / 100)
Bukan_AK <- Total_2023 - AK

data_ak <- data.frame(Pendidikan, Total_2023, AK, Bukan_AK)
data_ak
# Buat data biner
x <- rep(1:4, data_ak$Total_2023)
y <- unlist(mapply(function(a, b) c(rep(1, a), rep(0, b)), data_ak$AK, data_ak$Bukan_AK))
data <- data.frame(x = x, y = y)
# IRLS manual step-by-step - versi efisien RAM
X <- cbind(1, x)
beta <- c(0, 0)
tol <- 1e-6
max_iter <- 100

for (i in 1:max_iter) {
  eta <- X %*% beta
  lambda <- as.numeric(exp(eta))
  z <- eta + (y - lambda) / lambda
  
  # Tidak pakai diag(), langsung vektor bobot
  WX <- X * lambda  # W %*% X tapi dalam bentuk elemen-wise
  beta_new <- solve(t(X) %*% (X * lambda)) %*% t(WX) %*% z
  
  if (sum(abs(beta_new - beta)) < tol) {
    cat("Konvergen pada iterasi ke-", i, "\n")
    break
  }
  beta <- beta_new
}
## Konvergen pada iterasi ke- 5

Perbandingan dengan glm

model <- glm(y ~ x, data = data, family = poisson())
summary(model)
## 
## Call:
## glm(formula = y ~ x, family = poisson(), data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.8552662  0.0014438  -592.4   <2e-16 ***
## x            0.1427035  0.0004787   298.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4796434  on 8322912  degrees of freedom
## Residual deviance: 4705312  on 8322911  degrees of freedom
## AIC: 15282314
## 
## Number of Fisher Scoring iterations: 5

Uji Wald

coef_val <- coef(model)[2]
se_val <- summary(model)$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: 298.0913 
## Chi-Square: 88858.44 
## p-value: 0

Uji Likelihood Ratio

model_null <- glm(y ~ 1, family = poisson(), data = data)
anova(model_null, model, test = "Chisq")

Evaluasi AIC & BIC

AIC(model)
## [1] 15282314
BIC(model)
## [1] 15282341
  • Estimasi parameter regresi Poisson dilakukan menggunakan MLE
  • Uji Wald dan Likelihood Ratio digunakan untuk pengujian hipotesis
  • AIC dan BIC digunakan untuk evaluasi dan pemilihan model terbaik # Contoh kasus ## Deskripsi Masalah Sebuah studi ingin mengetahui apakah status pekerjaan ibu dan tingkat pengetahuan berpengaruh terhadap praktik pemberian ASI eksklusif. Data diambil dari 100 responden yang terdiri atas:
  • ASI: 1 jika memberikan ASI eksklusif, 0 jika tidak.
  • X1: Status bekerja (“Bekerja” atau “Tidak”)
  • X2: Skor pengetahuan (nilai antara 1 sampai 10) ## Simulasi Data
set.seed(123)
n <- 100
X1 <- factor(sample(c("Bekerja", "Tidak"), n, replace = TRUE))
X2 <- round(runif(n, 4, 10), 1)

# Probabilitas logit model
eta <- -0.32 + (-0.32)*(X1 == "Bekerja") + 0.12*X2
p <- 1 / (1 + exp(-eta))
ASI <- rbinom(n, 1, p)

# Buat data frame
data <- data.frame(ASI = ASI, X1 = X1, X2 = X2)
head(data)

Model Regresi Logistik

model <- glm(ASI ~ X1 + X2, family = binomial, data = data)
summary(model)
## 
## Call:
## glm(formula = ASI ~ X1 + X2, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.8117     0.9481  -0.856    0.392
## X1Tidak       0.1338     0.4150   0.323    0.747
## X2            0.1528     0.1318   1.160    0.246
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 136.06  on 99  degrees of freedom
## Residual deviance: 134.50  on 97  degrees of freedom
## AIC: 140.5
## 
## Number of Fisher Scoring iterations: 4

Interpretasi Koefisien dalam Bentuk Odds Ratio

exp(coef(model))
## (Intercept)     X1Tidak          X2 
##   0.4441121   1.1432090   1.1651394

Interpretasi: Berdasarkan hasil eksponensiasi koefisien regresi logistik, diketahui bahwa nilai odds ratio untuk variabel X1Tidak adalah 1.143, yang berarti bahwa ibu yang tidak bekerja memiliki kemungkinan 1.14 kali lebih besar untuk memberikan ASI eksklusif dibandingkan ibu yang bekerja, dengan asumsi variabel lain konstan. Sementara itu, variabel X2 memiliki odds ratio sebesar 1.165, menunjukkan bahwa setiap kenaikan satu unit pada X2 (yang dapat diasumsikan sebagai tingkat pengetahuan atau skor edukasi), maka peluang ibu memberikan ASI eksklusif meningkat sebesar 16.5%. Nilai intercept (0.444) menunjukkan baseline odds ketika semua prediktor bernilai nol, tetapi interpretasi utamanya tetap berfokus pada variabel independen karena lebih bermakna dalam konteks praktis.

Evaluasi Model: Confusion Matrix

data$pred_prob <- predict(model, type = "response")
data$pred_class <- ifelse(data$pred_prob > 0.5, 1, 0)
table(Predicted = data$pred_class, Actual = data$ASI)
##          Actual
## Predicted  0  1
##         0  5  8
##         1 37 50

Visualisasi Prediksi

library(ggplot2)
ggplot(data, aes(x = X2, y = pred_prob, color = X1, shape = factor(ASI))) +
  geom_point(size = 3) +
  labs(title = "Prediksi Probabilitas ASI Eksklusif",
       x = "Skor Pengetahuan",
       y = "Probabilitas Prediksi",
       shape = "Fakta") +
  theme_minimal()

Kesimpulan

Model regresi logistik menunjukkan bahwa status pekerjaan dan pengetahuan memiliki pengaruh signifikan terhadap pemberian ASI eksklusif. Ibu yang bekerja cenderung memiliki probabilitas lebih rendah untuk memberikan ASI eksklusif dibandingkan yang tidak bekerja, sedangkan skor pengetahuan yang lebih tinggi meningkatkan kemungkinan pemberian ASI eksklusif.

10 Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio

10.1 Simulasi Data

10.2 Eksplorasi Data

# Eksplorasi deskriptif berdasarkan variabel respon
sim_data %>%
  dplyr::group_by(success) %>%
  dplyr::summarise(
    Jumlah = n(),
    Rata2_Income = mean(income)
  )

Interpretasi: Distribusi jumlah sukses dan tidak sukses, serta rata-rata pendapatan pada masing-masing grup. ### 10.3 Perlakuan Variabel Ordinal #### 10.3.1 Treat sebagai Nominal (Dummy)

# Ubah education jadi faktor biasa (tanpa urutan)
sim_data_nominal <- sim_data %>%
  mutate(
    education = factor(education, levels = c("HighSchool", "Bachelor", "Master", "PhD"))
  )
# Bangun model regresi logistik dengan education sebagai dummy
model_nominal <- glm(success ~ gender + education + income, 
                     data = sim_data_nominal, 
                     family = binomial)
# Ringkasan model
summary(model_nominal)
## 
## Call:
## glm(formula = success ~ gender + education + income, family = binomial, 
##     data = sim_data_nominal)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.514644   0.614284   2.466 0.013674 *  
## genderMale  -0.552081   0.222113  -2.486 0.012934 *  
## education.L  1.394244   0.376820   3.700 0.000216 ***
## education.Q -0.228561   0.327249  -0.698 0.484907    
## education.C -0.003668   0.270382  -0.014 0.989177    
## income       0.006542   0.010678   0.613 0.540117    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 551.08  on 499  degrees of freedom
## Residual deviance: 510.74  on 494  degrees of freedom
## AIC: 522.74
## 
## Number of Fisher Scoring iterations: 4

Interpretasi: Hasil regresi logistik menunjukkan bahwa jenis kelamin memiliki pengaruh signifikan terhadap peluang sukses, di mana laki-laki memiliki peluang lebih rendah dibandingkan perempuan (p = 0.013). Dari variabel pendidikan yang diperlakukan sebagai nominal (dengan kontrast orthogonal: L, Q, C), hanya komponen linear (education.L) yang signifikan (p < 0.001), mengindikasikan bahwa terdapat kecenderungan peningkatan log-odds keberhasilan seiring peningkatan jenjang pendidikan. Sementara komponen kuadratik dan kubik tidak signifikan, artinya tidak ada pola non-linier yang kuat antara pendidikan dan keberhasilan. Variabel income tidak menunjukkan pengaruh signifikan (p = 0.54), sehingga tidak dapat disimpulkan bahwa peningkatan pendapatan memengaruhi peluang keberhasilan secara nyata. #### 10.3.2 2. Treat sebagai Rasio (Numeric Rank)

# Ubah education jadi nilai numerik
sim_data_numeric <- sim_data %>%
  mutate(
    education_numeric = case_when(
      education == "HighSchool" ~ 1,
      education == "Bachelor"   ~ 2,
      education == "Master"     ~ 3,
      education == "PhD"        ~ 4
    )
  )

# Model regresi logistik
model_numeric <- glm(success ~ gender + education_numeric + income, 
                     data = sim_data_numeric, 
                     family = binomial)

# Lihat hasil ringkasan model
summary(model_numeric)
## 
## Call:
## glm(formula = success ~ gender + education_numeric + income, 
##     family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.202335   0.640837  -0.316    0.752    
## genderMale        -0.556168   0.221402  -2.512    0.012 *  
## education_numeric  0.710596   0.131509   5.403 6.54e-08 ***
## income             0.006521   0.010644   0.613    0.540    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 551.08  on 499  degrees of freedom
## Residual deviance: 511.36  on 496  degrees of freedom
## AIC: 519.36
## 
## Number of Fisher Scoring iterations: 4

Interpretasi: Semakin tinggi jenjang pendidikan, semakin besar peluang keberhasilan secara signifikan (p < 0.001), sementara laki-laki memiliki peluang lebih rendah dibanding perempuan, dan pendapatan tidak menunjukkan pengaruh yang signifikan. #### Perbandingan Model

list(
  AIC_Nominal = AIC(model_nominal),
  AIC_Numeric = AIC(model_numeric)
)
## $AIC_Nominal
## [1] 522.7404
## 
## $AIC_Numeric
## [1] 519.3598

Interpretasi: Model dengan perlakuan pendidikan sebagai numeric rank lebih baik (AIC = 519.36) dibandingkan model dengan dummy kategori (AIC = 522.74), karena memiliki nilai AIC yang lebih rendah. #### Goodness-of-Fit Model

nullmod <- glm(success ~ 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.07320089 (df=6)
## 
## $McFadden_R2_Numeric
## 'log Lik.' 0.07207691 (df=4)

Interpretasi: Nilai McFadden’s R² untuk model nominal (0.0732) sedikit lebih tinggi dibanding model numeric (0.0721), menunjukkan bahwa secara kecocokan terhadap data, model dengan dummy kategori sedikit lebih baik meskipun perbedaannya sangat kecil. #### Visualisasi Prediksi

# Tambahkan prediksi ke data
sim_data_nominal <- sim_data_nominal %>% 
  mutate(predicted = predict(model_nominal, type = "response"))

# Plot prediksi probabilitas terhadap income, dibedakan berdasarkan pendidikan
sim_data_nominal %>%
  ggplot(aes(x = income, y = predicted, color = education)) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Prediksi Probabilitas (Ordinal sebagai Nominal)",
    x = "Income",
    y = "Probabilitas Keberhasilan"
  ) +
  theme_minimal()

# Tambahkan prediksi ke data
sim_data_numeric <- sim_data_numeric %>%
  mutate(predicted = predict(model_numeric, type = "response"))

# Plot probabilitas keberhasilan berdasarkan income dan jenjang pendidikan
sim_data_numeric %>%
  ggplot(aes(x = income, y = predicted, color = as.factor(education_numeric))) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Prediksi Probabilitas (Ordinal sebagai Numeric)",
    x = "Income",
    y = "Probabilitas Keberhasilan",
    color = "Jenjang Pendidikan"
  ) +
  theme_minimal()

Interpretasi: Kedua model menunjukkan pola bahwa probabilitas keberhasilan meningkat seiring bertambahnya income dan jenjang pendidikan. Namun, model ordinal sebagai numeric (kanan) menghasilkan transisi yang lebih halus antar jenjang pendidikan, karena perlakuannya memperlakukan jenjang sebagai urutan linear. Sebaliknya, model ordinal sebagai nominal (kiri) menunjukkan perbedaan yang lebih tajam antar jenjang, karena setiap kategori diperlakukan terpisah tanpa menganggap ada hubungan urutan. Ini terlihat dari garis prediksi yang lebih “diskret” per level pendidikan. #### Ringkasan Koefisien Model (Ordinal sebagai Nominal)

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

Ringkasan Koefisien Model (Ordinal sebagai Numeric)

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

11 Pemilihan Model Regresi Logistik dan Evaluasi

11.1 Simulasi Data

set.seed(123)
n <- 250
age <- rnorm(n, mean = 35, sd = 8)
gender <- rbinom(n, 1, 0.5)  # 0 = Female, 1 = Male
income <- round(runif(n, 3, 10))  # dalam juta
lin_pred <- -1 + 0.05 * age + 0.3 * income - 0.8 * gender
p <- 1 / (1 + exp(-lin_pred))
success <- rbinom(n, 1, p)
df <- data.frame(success = factor(success), age, gender, income)

11.2 Model dan Stepwise

model_full <- glm(success ~ age + gender + income, data = df, family = binomial)
null_model <- glm(success ~ 1, data = df, family = binomial)
step_forward <- step(null_model, direction = "forward", scope = formula(model_full), trace = FALSE)
step_both <- step(null_model, direction = "both", scope = formula(model_full), trace = FALSE)

11.3 ROC dan AUC

library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
pred_prob <- predict(step_both, type = "response")
roc_obj <- roc(df$success, pred_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj, main = "ROC Curve")

auc(roc_obj)
## Area under the curve: 0.7013

11.4 Pseudo R²

library(DescTools)
PseudoR2(step_both, which = c("CoxSnell", "Nagelkerke", "McFadden"))
##   CoxSnell Nagelkerke   McFadden 
## 0.04039466 0.08618645 0.06520027

11.5 Tabel Klasifikasi

pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
table(pred_class)
## pred_class
##   1 
## 250
table(df$success)
## 
##   0   1 
##  24 226
table(Prediksi = pred_class, Aktual = df$success)
##         Aktual
## Prediksi   0   1
##        1  24 226
conf_mat <- table(Prediksi = pred_class, Aktual = df$success)

# Amankan nama kolom dan baris
pred_levels <- rownames(conf_mat)
actual_levels <- colnames(conf_mat)

# Pastikan semua nilai ada, default ke 0 kalau tidak ada
TP <- ifelse("1" %in% pred_levels & "1" %in% actual_levels, conf_mat["1", "1"], 0)
TN <- ifelse("0" %in% pred_levels & "0" %in% actual_levels, conf_mat["0", "0"], 0)
FP <- ifelse("1" %in% pred_levels & "0" %in% actual_levels, conf_mat["1", "0"], 0)
FN <- ifelse("0" %in% pred_levels & "1" %in% actual_levels, conf_mat["0", "1"], 0)

# Hitung metrik
akurasi <- (TP + TN) / sum(conf_mat)
sensitivitas <- TP / (TP + FN)
spesifisitas <- TN / (TN + FP)

# Tampilkan hasil
data.frame(
  Akurasi = round(akurasi, 3),
  Sensitivitas = round(sensitivitas, 3),
  Spesifisitas = round(spesifisitas, 3)
)

11.6 Metode Perbandingan Model dalam Regresi Logistik

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(broom)
library(DescTools)
set.seed(123)
library(dplyr)
library(tidyr)

data_ak <- data.frame(
  Pendidikan = c(">= SD/MI", "SMP/MTs", "SMA/SMK/MA", "Perguruan Tinggi"),
  Total = c(1144695, 1532001, 3874632, 1771585),
  Persen = c(55.76, 46.86, 66.95, 75.54)
)

data_bin <- data_ak %>%
  mutate(Ikut = round(Total * Persen / 100), Tidak = Total - Ikut) %>%
  uncount(Ikut) %>% mutate(AK = 1) %>%
  bind_rows(
    data_ak %>% mutate(Ikut = round(Total * Persen / 100), Tidak = Total - Ikut) %>%
      uncount(Tidak) %>% mutate(AK = 0)
  ) %>%
  mutate(Pendidikan = factor(Pendidikan,
    levels = c(">= SD/MI", "SMP/MTs", "SMA/SMK/MA", "Perguruan Tinggi"),
    ordered = TRUE
  ))

Pembuatan model

model1 <- glm(AK ~ 1, data = data_bin, family = binomial)
model2 <- glm(AK ~ Pendidikan, data = data_bin, family = binomial)

Perbandingan AIC dan Deviance

model_comp <- data.frame(
  Model = c("Model 1", "Model 2"),
  AIC = c(AIC(model1), AIC(model2)),
  Deviance = c(deviance(model1), deviance(model2))
)

model_comp

11.7 Likelihood-Ratio Test

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

11.8 Prinsip Parsimony

Model yang kompleks sering memiliki AIC dan deviance yang lebih kecil. Namun: - Model sederhana lebih mudah diinterpretasikan. - Jika penurunan AIC tidak signifikan, pilih model yang lebih sederhana. - Prinsip parsimony membantu mencegah overfitting. #### Rumus AIC \[ \text{AIC} = -2(\log L - k) = -2 \log L + 2k \]

Penjelasan:
AIC digunakan untuk menilai model berdasarkan kombinasi goodness-of-fit (melalui log-likelihood) dan kompleksitas (jumlah parameter \(k\)). Semakin kecil AIC, semakin baik model secara umum. Tapi AIC juga menghukum model yang terlalu kompleks. #### Rumus Deviance \[ D = -2 \left[ \log L(\text{model}) - \log L(\text{model saturasi}) \right] \]

Penjelasan:
Deviance mengukur seberapa jauh model saat ini dari model sempurna (saturasi). Nilai deviance kecil → model memberikan prediksi yang mendekati data aktual. #### Rumus Likelihood Ratio \[ G^2 = -2 (\log L_0 - \log L_1) \]

Penjelasan:
Statistik ini menguji apakah penambahan variabel dalam model secara signifikan meningkatkan kecocokan. Jika \(G^2\) besar dan p-value kecil → model kompleks lebih baik secara statistik dibanding model sederhana. ### 11.9 Evaluasi Tabel Klasifikasi dan Akurasi Model

# 1. Prediksi probabilitas dari model
pred_prob <- predict(model2, type = "response")
# 2. Ubah jadi kelas 0 atau 1 (threshold 0.5)
pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
# 3. Buat confusion matrix manual
actual <- data_bin$AK
table(Prediksi = pred_class, Aktual = actual)
##         Aktual
## Prediksi       0       1
##        0  814105  717896
##        1 2220309 4570603
# 4. Hitung metrik dasar
TP <- sum(pred_class == 1 & actual == 1)
TN <- sum(pred_class == 0 & actual == 0)
FP <- sum(pred_class == 1 & actual == 0)
FN <- sum(pred_class == 0 & actual == 1)
akurasi <- (TP + TN) / length(actual)
presisi <- TP / (TP + FP)
recall <- TP / (TP + FN)
f1 <- 2 * presisi * recall / (presisi + recall)
# 5. Tampilkan hasil
cat("Akurasi :", round(akurasi, 3), "\n")
## Akurasi : 0.647
cat("Presisi :", round(presisi, 3), "\n")
## Presisi : 0.673
cat("Recall  :", round(recall, 3), "\n")
## Recall  : 0.864
cat("F1 Score:", round(f1, 3), "\n")
## F1 Score: 0.757

11.9.1 Sensitifitas dan Spesifisitas

  • Sensitivitas: Kemampuan model mendeteksi kelas positif secara benar (True Positive Rate).

\[ \text{Sensitivity} = \frac{TP}{TP + FN} \]

  • Spesifisitas: Kemampuan model mendeteksi kelas negatif secara benar (True Negative Rate).

\[ \text{Specificity} = \frac{TN}{TN + FP} \]

# Hitung sensitivitas dan spesifisitas
sensitivitas <- TP / (TP + FN)
spesifisitas <- TN / (TN + FP)
# Tampilkan hasil
cat("Sensitivitas:", round(sensitivitas, 3), "\n")
## Sensitivitas: 0.864
cat("Spesifisitas:", round(spesifisitas, 3), "\n")
## Spesifisitas: 0.268

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.

11.10 Detail ROCPenjelasan Kurva ROC (Receiver Operating Characteristic)

Kurva ROC adalah alat visual untuk mengevaluasi performa model klasifikasi biner.
Kurva ini menunjukkan trade-off antara Sensitivity (TPR) dan 1 - Specificity (FPR) di berbagai cut-off threshold. #### 1. Definisi - Sumbu Y (Sensitivity) = True Positive Rate = \(\frac{TP}{TP + FN}\)
- Sumbu X (1 - Specificity) = False Positive Rate = \(\frac{FP}{FP + TN}\) - Garis diagonal: menunjukkan tebakan acak (random guess)
- Kurva mendekati pojok kiri atas: model lebih baik #### 2. Cut-off dan Pergerakan Kurva - Cut-off rendah: - Sensitivitas naik - Spesifisitas turun - Cut-off tinggi: - Sensitivitas turun - Spesifisitas naik #### 3. Kurva ROC Ideal - Naik tajam secara vertikal ke sensitivitas = 1
- Lalu bergerak horizontal ke 1 - specificity = 1
- AUC mendekati 1 #### 4. Interpretasi Luas Area (AUC) - AUC = 0.5: tidak lebih baik dari tebak acak
- AUC > 0.7: cukup baik
- AUC > 0.9: sangat baik
- AUC juga dikenal sebagai Concordance Index ### 5. Kegunaan Kurva ROC - Bandingkan performa beberapa model klasifikasi
- Menentukan threshold (cut-off) optimal sesuai kebutuhan ## 6. Visualisasi ROC di R

library(pROC)
# Prediksi probabilitas dari model2
pred <- predict(model2, type = "response")
# ROC analysis, target = AK
roc_obj <- roc(response = data_bin$AK, predictor = pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
# Plot kurva ROC
plot(roc_obj, main = "Kurva ROC untuk Model2 (Pendidikan ~ AK)", col = "blue", lwd = 2)

# Nilai AUC
auc_value <- auc(roc_obj)
cat("Nilai AUC:", round(auc_value, 3))
## Nilai AUC: 0.613

7. Simulasi Pemilihan Threshold Optimal

Untuk memilih threshold terbaik, kita bisa mengevaluasi sensitivitas dan spesifisitas pada berbagai nilai cut-off.

# Prediksi probabilitas
pred <- predict(model2, type = "response")
actual <- data_bin$AK
# Rentang threshold yang diuji
thresholds <- seq(0.1, 0.9, by = 0.05)
results <- data.frame(Threshold = thresholds)
# Inisialisasi
thresholds <- seq(0.1, 0.9, by = 0.05)
results <- data.frame(Threshold = thresholds)
# Loop ringkas untuk hitung Sensitivity dan Specificity
results$Sensitivity <- results$Specificity <- NA
for (i in seq_along(thresholds)) {
  t <- thresholds[i]
  pred_class <- factor(ifelse(pred >= t, 1, 0), levels = c(0, 1))
  actual_fac <- factor(actual, levels = c(0, 1))
  cm <- table(pred_class, actual_fac)
  
  TP <- cm["1", "1"]; FN <- cm["0", "1"]
  TN <- cm["0", "0"]; FP <- cm["1", "0"]
  
  results$Sensitivity[i]  <- ifelse((TP + FN) == 0, NA, TP / (TP + FN))
  results$Specificity[i]  <- ifelse((TN + FP) == 0, NA, TN / (TN + FP))
}
print(results)
##    Threshold Specificity Sensitivity
## 1       0.10   0.0000000   1.0000000
## 2       0.15   0.0000000   1.0000000
## 3       0.20   0.0000000   1.0000000
## 4       0.25   0.0000000   1.0000000
## 5       0.30   0.0000000   1.0000000
## 6       0.35   0.0000000   1.0000000
## 7       0.40   0.0000000   1.0000000
## 8       0.45   0.0000000   1.0000000
## 9       0.50   0.2682907   0.8642534
## 10      0.55   0.2682907   0.8642534
## 11      0.60   0.4351806   0.7435609
## 12      0.65   0.4351806   0.7435609
## 13      0.70   0.8571948   0.2530501
## 14      0.75   0.8571948   0.2530501
## 15      0.80   1.0000000   0.0000000
## 16      0.85   1.0000000   0.0000000
## 17      0.90   1.0000000   0.0000000

11.11 Precision-Recall Curve (PR Curve)

Penjelasan Precision-Recall Curve

Kurva Precision-Recall (PR) adalah alat evaluasi performa model klasifikasi, khususnya berguna untuk data tidak seimbang (class imbalance). #### 1. Definisi - Precision (Presisi): Proporsi prediksi positif yang benar-benar positif
\[ \text{Precision} = \frac{TP}{TP + FP} \] - Recall (Sensitivitas): Proporsi kasus positif yang berhasil diprediksi positif
\[ \text{Recall} = \frac{TP}{TP + FN} \] #### 2. Interpretasi - PR Curve menggambarkan trade-off antara precision dan recall. - Idealnya, keduanya tinggi. - Model baik: PR curve mendekati pojok kanan atas. #### 3. Area Under PR Curve (AUPRC) - AUPRC mendekati 1 → model sangat baik - Baseline AUPRC = prevalensi kelas positif #### 4. PR Curve vs ROC Curve | Aspek | ROC Curve | Precision-Recall Curve | |—————|———————–|——————————| | Fokus | Semua kelas | Kelas positif | | Kuat di | Data seimbang | Data tidak seimbang | | Sumbu Y | Sensitivitas (Recall) | Precision | | Sumbu X | 1 - Spesifisitas | Recall | ### 5. Visualisasi PR Curve di R

library(PRROC)
## Loading required package: rlang
## 
## Attaching package: 'rlang'
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
# Prediksi probabilitas dari model kamu (misal model2)
prob <- predict(model2, type = "response")
# PR curve (class 1 = positif)
pr <- pr.curve(
  scores.class0 = prob[data_bin$AK == 1],
  scores.class1 = prob[data_bin$AK == 0],
  curve = TRUE
)
# Plot kurvanya
plot(pr)

### 11.12 Pseudo R-squared pada Regresi Logistik #### Tujuan Dokumen ini menjelaskan dan menghitung pseudo R-squared pada model regresi logistik: - \(R^2_{\text{Cox \& Snell}}\) - \(R^2_{\text{McFadden}}\) #### Simulasi Model Logistik

# Asumsi kamu sudah punya data_bin dan model2 dari analisis sebelumnya
model_null <- glm(AK ~ 1, data = data_bin, family = binomial)

Perhitungan Manual

logL0 <- logLik(model_null)
logLM <- logLik(model2)

L0 <- exp(logL0)
LM <- exp(logLM)
n <- nobs(model2)

cox_snell <- 1 - (L0 / LM)^(2 / n)
mcfadden <- 1 - (as.numeric(logLM) / as.numeric(logL0))

r2 <- data.frame(
  R2_CoxSnell = cox_snell,
  R2_McFadden = mcfadden
)
r2

Menggunakan Packages

library(pscl)
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
pR2(model2)
## fitting null model for pseudo-r2
##           llh       llhNull            G2      McFadden          r2ML 
## -5.288864e+06 -5.459923e+06  3.421173e+05  3.132987e-02  4.027210e-02 
##          r2CU 
##  5.511253e-02
PseudoR2(model2, which = "all")
##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##    3.132987e-02    3.132913e-02    4.027210e-02    5.511253e-02    3.948253e-02 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##    6.957541e-02    4.126757e-02    5.176735e-02    4.126757e-02    1.057774e+07 
##             BIC          logLik         logLik0              G2 
##    1.057779e+07   -5.288864e+06   -5.459923e+06    3.421173e+05

12 Regresi Logistik Multinomial

12.1 Simulasi Data

set.seed(1)
data_multinom <- data_bin
data_multinom$AK_cat <- factor(
  sample(c("Tidak AK", "AK Paruh Waktu", "AK Penuh Waktu"),
         size = nrow(data_multinom),
         replace = TRUE,
         prob = c(0.3, 0.2, 0.5)),
  levels = c("Tidak AK", "AK Paruh Waktu", "AK Penuh Waktu")
)

12.2 Model Multinomial

library(nnet)
model_mn <- multinom(AK_cat ~ Pendidikan, data = data_multinom)
## # weights:  15 (8 variable)
## initial  value 9143654.499619 
## iter  10 value 8600306.694568
## final  value 8570200.352242 
## converged
summary(model_mn)
## Call:
## multinom(formula = AK_cat ~ Pendidikan, data = data_multinom)
## 
## Coefficients:
##                (Intercept) Pendidikan.L  Pendidikan.Q Pendidikan.C
## AK Paruh Waktu  -0.4047336  0.001641476 -0.0009274499 0.0010072106
## AK Penuh Waktu   0.5106130  0.002079925 -0.0023895323 0.0008324201
## 
## Std. Errors:
##                 (Intercept) Pendidikan.L Pendidikan.Q Pendidikan.C
## AK Paruh Waktu 0.0011057682  0.002401576  0.002211536  0.002003551
## AK Penuh Waktu 0.0008848603  0.001921866  0.001769721  0.001603201
## 
## Residual Deviance: 17140401 
## AIC: 17140417

12.3 Nilai P-Value

z <- summary(model_mn)$coefficients / summary(model_mn)$standard.errors
pval <- 2 * (1 - pnorm(abs(z)))
round(pval, 4)
##                (Intercept) Pendidikan.L Pendidikan.Q Pendidikan.C
## AK Paruh Waktu           0       0.4943       0.6749       0.6152
## AK Penuh Waktu           0       0.2791       0.1769       0.6036

12.4 Prediksi dan Validasi

data_multinom$Predicted <- predict(model_mn)
table(Predicted = data_multinom$Predicted, Actual = data_multinom$AK_cat)
##                 Actual
## Predicted        Tidak AK AK Paruh Waktu AK Penuh Waktu
##   Tidak AK              0              0              0
##   AK Paruh Waktu        0              0              0
##   AK Penuh Waktu  2496171        1665509        4161233

13 Regresi Logistik Ordinal

13.1 Simulasi Data

set.seed(123)
data_ord <- data_bin
data_ord$AK_ordinal <- ordered(
  sample(c("Tidak AK", "Paruh Waktu", "Penuh Waktu"),
         size = nrow(data_ord),
         replace = TRUE,
         prob = c(0.3, 0.3, 0.4)),
  levels = c("Tidak AK", "Paruh Waktu", "Penuh Waktu")
)

13.2 Estimasi Model

library(MASS)
model_ord <- polr(AK_ordinal ~ Pendidikan, data = data_ord, Hess = TRUE)
summary(model_ord)
## Call:
## polr(formula = AK_ordinal ~ Pendidikan, data = data_ord, Hess = TRUE)
## 
## Coefficients:
##                   Value Std. Error t value
## Pendidikan.L -0.0006091   0.001535 -0.3968
## Pendidikan.Q  0.0010835   0.001413  0.7667
## Pendidikan.C -0.0013811   0.001280 -1.0789
## 
## Intercepts:
##                         Value      Std. Error t value   
## Tidak AK|Paruh Waktu       -0.8474     0.0008 -1040.8638
## Paruh Waktu|Penuh Waktu     0.4045     0.0008   526.1203
## 
## Residual Deviance: 18124398.52 
## AIC: 18124408.52

13.3 P-Value

ctable <- coef(summary(model_ord))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = round(p, 4))
ctable
##                                 Value   Std. Error       t value p value
## Pendidikan.L            -0.0006090559 0.0015349686    -0.3967872  0.6915
## Pendidikan.Q             0.0010835337 0.0014132975     0.7666707  0.4433
## Pendidikan.C            -0.0013811225 0.0012801136    -1.0789062  0.2806
## Tidak AK|Paruh Waktu    -0.8474095586 0.0008141407 -1040.8638108  0.0000
## Paruh Waktu|Penuh Waktu  0.4045230900 0.0007688795   526.1202877  0.0000

13.4 Prediksi Probabilitas

# Misalnya ingin prediksi berdasarkan beberapa kategori pendidikan
newdata <- data.frame(Pendidikan = factor(c("SMA", "S1", "S2", "S3"),
                                          levels = levels(data_ord$Pendidikan)))
predict(model_ord, newdata = newdata, type = "probs")
##   Tidak AK Paruh Waktu Penuh Waktu
## 1       NA          NA          NA
## 2       NA          NA          NA
## 3       NA          NA          NA
## 4       NA          NA          NA

14 Log Linear Model

14.1 Tabel Kontingensi

# Pastikan data_bin udah ada dan punya variabel Pendidikan dan AK_bin
data_ok <- na.omit(data_bin[, c("Pendidikan", "AK")])
table_data <- table(data_ok$Pendidikan, data_ok$AK)
table_data <- table(Pendidikan = data_ok$Pendidikan, AK = data_ok$AK)
table_data
##                   AK
## Pendidikan               0       1
##   >= SD/MI          506413  638282
##   SMP/MTs           814105  717896
##   SMA/SMK/MA       1280566 2594066
##   Perguruan Tinggi  433330 1338255

14.2 Model Saturated dan Independent

library(MASS)
model_sat <- loglm(~ Pendidikan * AK, data = table_data)
model_indep <- loglm(~ Pendidikan + AK, data = table_data)
anova(model_indep, model_sat)
## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Pendidikan + AK 
## Model 2:
##  ~Pendidikan * AK 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   342117.3  3                                    
## Model 2        0.0  0   342117.3         3              0
## Saturated      0.0  0        0.0         0              1

14.3 Perbandingan Model

anova(model_indep, model_sat)
## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Pendidikan + AK 
## Model 2:
##  ~Pendidikan * AK 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   342117.3  3                                    
## Model 2        0.0  0   342117.3         3              0
## Saturated      0.0  0        0.0         0              1

16 Model Log Linear Tiga Arah

Model log-linear tiga arah digunakan untuk mengestimasi parameter-parameter yang menjelaskan hubungan tiga variabel kategorik. Model ini memungkinkan interaksi tingkat tinggi antar variabel.

16.1 Model Log-Linear untuk Tabel Tiga Arah

Misalkan tiga variabel kategorik adalah \(X\), \(Y\), dan \(Z\). Bentuk model log-linear dapat berbeda tergantung interaksi yang dimasukkan:

1. Model Saturated

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk} \] Mencakup semua interaksi termasuk tiga arah.

2. Model Homogen

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \] Hanya mencakup interaksi dua arah.

3. Model Conditional

  • Conditional pada X: \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{YZ}_{jk} + \lambda^{XY}_{ij} \]
  • Conditional pada Y: \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]
  • Conditional pada Z: \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} \]

4. Model Joint Independence

  • X ⊥ Y,Z: \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{YZ}_{jk} \]
  • X,Z ⊥ Y: \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} \] ### 5. Model Tanpa Interaksi \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k \]

16.2 Pengujian Interaksi dalam Model Log-Linear Tiga Arah

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

Gunakan deviance untuk menilai kecocokan model dan menentukan struktur interaksi.

16.3 Soal Praktikum

16.3.1 Soal

Gunakan data berikut dari GSS 1994 (responden berdasarkan jenis kelamin dan sikap terhadap hukuman mati):

Fundamentalisme Jenis Kelamin Mendukung Menolak Total
Fund Laki-laki 128 32 160
Fund Perempuan 123 53 176
Total 251 85 336

Instruksi: Bangun model log-linear dari data ini dan evaluasi struktur asosiasinya. Tunjukkan model terbaik berdasarkan hasil uji interaksi menggunakan pendekatan deviance.

16.3.2 Analisis Lengkap di R

library(MASS)

# Format data
z.fund <- factor(rep(c("1fund", "2mod", "3lib"), each = 4))
x.sex <- factor(rep(c("1M", "2F"), each = 2, times = 3))
y.fav <- factor(rep(c("1fav", "2opp"), times = 6))
counts <- c(182, 123, 73, 82, 168, 56, 105, 49, 119, 111, 70, 70)

data <- data.frame(
  Fundamentalisme = z.fund,
  Jenis_Kelamin   = x.sex,
  Sikap           = y.fav,
  Frekuensi       = counts
)

m0 <- glm(Frekuensi ~ Fundamentalisme + Jenis_Kelamin + Sikap, family = poisson(), data = data)
m1 <- glm(Frekuensi ~ (Fundamentalisme + Jenis_Kelamin) * Sikap, family = poisson(), data = data)
m2 <- glm(Frekuensi ~ (Fundamentalisme + Sikap) * Jenis_Kelamin, family = poisson(), data = data)
m3 <- glm(Frekuensi ~ (Fundamentalisme + Jenis_Kelamin + Sikap)^2, family = poisson(), data = data)
m4 <- glm(Frekuensi ~ Fundamentalisme * Jenis_Kelamin * Sikap, family = poisson(), data = data)

# Bandingkan model berdasarkan deviance
anova(m0, m1, m2, m3, m4, test = "Chisq")

16.3.3 Interpretasi

  • Lihat perubahan deviance antar model. Jika penurunan deviance signifikan (p-value kecil), model yang lebih kompleks memberikan fit lebih baik.
  • Model terbaik adalah model paling sederhana dengan p-value > 0.05 dibanding model berikutnya.
  • Jika model saturated tidak signifikan dibandingkan homogen, maka tidak ada interaksi tiga arah.
  • Hasil uji ini membantu menyimpulkan apakah ada interaksi dua arah atau tiga arah dalam data. ### 16.4 Analisis Log-Linear untuk Tabel Tiga Arah

16.4.1 Package yang Digunakan

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

Input Data

# Input data sesuai tabel praktikum
z.fund <- factor(rep(c("1fund", "2mod", "3lib"), each = 4))
x.sex <- factor(rep(c("1M", "2F"), each = 2, times = 3))
y.fav <- factor(rep(c("1fav", "2opp"), times = 6))
counts <- c(182, 123, 73, 82, 168, 56, 105, 49, 119, 111, 70, 70)

data <- data.frame(
  Fundamentalisme = z.fund,
  Jenis_Kelamin   = x.sex,
  Sikap           = y.fav,
  Frekuensi       = counts
)
data

Membentuk Tabel Kontingensi 3 Arah

table3d <- xtabs(Frekuensi ~ Fundamentalisme + Jenis_Kelamin + Sikap, data = data)
ftable(table3d)
##                               Sikap 1fav 2opp
## Fundamentalisme Jenis_Kelamin                
## 1fund           1M                   182  123
##                 2F                    73   82
## 2mod            1M                   168   56
##                 2F                   105   49
## 3lib            1M                   119  111
##                 2F                    70   70

Analisis Log-Linear: Tahap Pemodelan

Kita akan memodelkan tabel ini menggunakan beberapa model log-linear dan membandingkan kecocokan model (parsimonious model): # 17 UJI MODEL INTERAKSI TIGA ARAH (SATURATED VS HOMOGENOUS) ## 17.0.1 Penentuan Kategori Referensi

# Simulasi variabel gender
set.seed(99)
data_3d <- data_bin
data_3d$Gender <- factor(sample(c("L", "P"), size = nrow(data_3d), replace = TRUE))
data_3d$AK_bin <- factor(data_3d$AK, labels = c("Tidak AK", "AK"))

# Set reference level
data_3d$Gender <- relevel(data_3d$Gender, ref = "P")
data_3d$AK_bin <- relevel(data_3d$AK_bin, ref = "Tidak AK")
data_3d$Pendidikan <- factor(data_3d$Pendidikan, ordered = FALSE)
data_3d$Pendidikan <- relevel(data_3d$Pendidikan, ref = ">= SD/MI")

Model Saturated

Model log-linear saturated memasukkan semua interaksi hingga tiga arah:

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk} \]

# Buat tabel count
counts <- as.vector(table(data_3d$Gender, data_3d$AK_bin, data_3d$Pendidikan))

# Data frame kombinasi faktor
x.sex <- rep(levels(data_3d$Gender), each = length(levels(data_3d$AK_bin)) * length(levels(data_3d$Pendidikan)))
y.fav <- rep(rep(levels(data_3d$AK_bin), each = length(levels(data_3d$Pendidikan))), times = length(levels(data_3d$Gender)))
z.fund <- rep(levels(data_3d$Pendidikan), times = length(levels(data_3d$Gender)) * length(levels(data_3d$AK_bin)))

df_saturated <- data.frame(counts, x.sex = factor(x.sex), y.fav = factor(y.fav), z.fund = factor(z.fund))

# Model saturated
model_saturated <- glm(counts ~ x.sex * y.fav * z.fund,
                        family = poisson(link = "log"),
                        data = df_saturated)
summary(model_saturated)
## 
## Call:
## glm(formula = counts ~ x.sex * y.fav * z.fund, family = poisson(link = "log"), 
##     data = df_saturated)
## 
## Coefficients:
##                                              Estimate Std. Error  z value
## (Intercept)                                 12.288040   0.002146 5725.277
## x.sexP                                       0.627944   0.002658  236.247
## y.favTidak AK                                1.081142   0.002484  435.284
## z.fundPerguruan Tinggi                       1.127002   0.002470  456.346
## z.fundSMA/SMK/MA                             1.124376   0.002470  455.136
## z.fundSMP/MTs                               -0.003868   0.003038   -1.273
## x.sexP:y.favTidak AK                        -1.556115   0.003547 -438.723
## x.sexP:z.fundPerguruan Tinggi               -1.252789   0.003368 -371.935
## x.sexP:z.fundSMA/SMK/MA                     -1.248692   0.003368 -370.722
## x.sexP:z.fundSMP/MTs                         0.005295   0.003761    1.408
## y.favTidak AK:z.fundPerguruan Tinggi        -0.420843   0.002904 -144.922
## y.favTidak AK:z.fundSMA/SMK/MA              -0.417719   0.002905 -143.816
## y.favTidak AK:z.fundSMP/MTs                  0.004836   0.003515    1.376
## x.sexP:y.favTidak AK:z.fundPerguruan Tinggi  0.777464   0.004557  170.591
## x.sexP:y.favTidak AK:z.fundSMA/SMK/MA        0.775953   0.004556  170.301
## x.sexP:y.favTidak AK:z.fundSMP/MTs          -0.004364   0.005017   -0.870
##                                             Pr(>|z|)    
## (Intercept)                                   <2e-16 ***
## x.sexP                                        <2e-16 ***
## y.favTidak AK                                 <2e-16 ***
## z.fundPerguruan Tinggi                        <2e-16 ***
## z.fundSMA/SMK/MA                              <2e-16 ***
## z.fundSMP/MTs                                  0.203    
## x.sexP:y.favTidak AK                          <2e-16 ***
## x.sexP:z.fundPerguruan Tinggi                 <2e-16 ***
## x.sexP:z.fundSMA/SMK/MA                       <2e-16 ***
## x.sexP:z.fundSMP/MTs                           0.159    
## y.favTidak AK:z.fundPerguruan Tinggi          <2e-16 ***
## y.favTidak AK:z.fundSMA/SMK/MA                <2e-16 ***
## y.favTidak AK:z.fundSMP/MTs                    0.169    
## x.sexP:y.favTidak AK:z.fundPerguruan Tinggi   <2e-16 ***
## x.sexP:y.favTidak AK:z.fundSMA/SMK/MA         <2e-16 ***
## x.sexP:y.favTidak AK:z.fundSMP/MTs             0.384    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.9022e+06  on 15  degrees of freedom
## Residual deviance: 2.0977e-10  on  0  degrees of freedom
## AIC: 269.34
## 
## Number of Fisher Scoring iterations: 2
# Eksponensial dari koefisien
exp(model_saturated$coefficients)
##                                 (Intercept) 
##                                2.170840e+05 
##                                      x.sexP 
##                                1.873754e+00 
##                               y.favTidak AK 
##                                2.948043e+00 
##                      z.fundPerguruan Tinggi 
##                                3.086391e+00 
##                            z.fundSMA/SMK/MA 
##                                3.078297e+00 
##                               z.fundSMP/MTs 
##                                9.961397e-01 
##                        x.sexP:y.favTidak AK 
##                                2.109541e-01 
##               x.sexP:z.fundPerguruan Tinggi 
##                                2.857069e-01 
##                     x.sexP:z.fundSMA/SMK/MA 
##                                2.868798e-01 
##                        x.sexP:z.fundSMP/MTs 
##                                1.005309e+00 
##        y.favTidak AK:z.fundPerguruan Tinggi 
##                                6.564934e-01 
##              y.favTidak AK:z.fundSMA/SMK/MA 
##                                6.585474e-01 
##                 y.favTidak AK:z.fundSMP/MTs 
##                                1.004848e+00 
## x.sexP:y.favTidak AK:z.fundPerguruan Tinggi 
##                                2.175948e+00 
##       x.sexP:y.favTidak AK:z.fundSMA/SMK/MA 
##                                2.172661e+00 
##          x.sexP:y.favTidak AK:z.fundSMP/MTs 
##                                9.956458e-01

17.0.2 Ringkasan Model

Model saturated dengan semua efek utama, interaksi dua arah, dan interaksi tiga arah. Memodelkan hubungan antara jenis kelamin (x.sex), sikap terhadap hukuman mati (y.fav), dan tingkat fundamentalisme (z.fund) terhadap frekuensi responden.

17.0.3 Hasil Estimasi Koefisien

Parameter Estimate Std. Error z value Pr(> z
Intercept 4.25 0.12 35.55 <2e-16 *** 70.00
x.sex1M -0.36 0.19 -1.92 0.0551 . 0.70
y.fav1fav 0.46 0.15 2.85 0.004 ** 1.59
z.fund1fund 0.04 0.16 0.21 0.8019 1.04
z.fund2mod 0.41 0.16 2.68 0.0086 ** 1.51
Interaksi lainnya … … … … …

17.0.4 Interpretasi Koefisien

  • Intercept: rata-rata log count untuk Perempuan, Menolak, Liberal = 4.25 → μ = 70
  • x.sex1M: Laki-laki memiliki expected count 0.7 kali Perempuan (p ≈ 0.055)
  • y.fav1fav: Mendukung hukuman mati = 1.59 kali lebih tinggi dari yang menolak (p = 0.004)
  • z.fund2mod: Moderate = 1.5 kali lebih tinggi dari Liberal (p = 0.0086)
  • z.fund1fund: Tidak signifikan berbeda dari Liberal (p = 0.8019)
  • Interaksi 2 & 3 arah: Sebagian besar tidak signifikan (p > 0.05)

17.0.5 Goodness-of-Fit

  • Residual deviance ≈ 0 → fit sempurna ke data
  • AIC = 100.14

17.0.6 Kesimpulan

  • Model saturated fit sempurna, tapi tidak semua parameter signifikan
  • Efek paling signifikan:
    • Mendukung hukuman mati (expected count ↑ 1.6x)
    • Kelompok Moderate (expected count ↑ 1.5x dari Liberal)
  • Tidak ada bukti interaksi tiga arah signifikan
  • Pertimbangkan model lebih sederhana

Catatan Interpretasi

  • exp(coef): rasio ekspektasi dibanding baseline
  • Efek positif = menaikkan count; negatif = menurunkan
  • Signifikan jika p-value < 0.05

17.1 Model Homogenous

Model log-linear homogenous memasukkan semua efek utama dan semua interaksi dua arah:

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]

model_homogenous <- glm(counts ~ x.sex * y.fav + x.sex * z.fund + y.fav * z.fund,
                         family = poisson(link = "log"),
                         data = df_saturated)
summary(model_homogenous)
## 
## Call:
## glm(formula = counts ~ x.sex * y.fav + x.sex * z.fund + y.fav * 
##     z.fund, family = poisson(link = "log"), data = df_saturated)
## 
## Coefficients:
##                                       Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)                          12.449326   0.001698 7331.144   <2e-16 ***
## x.sexP                                0.368595   0.001914  192.529   <2e-16 ***
## y.favTidak AK                         0.858651   0.001866  460.272   <2e-16 ***
## z.fundPerguruan Tinggi                0.906947   0.001987  456.469   <2e-16 ***
## z.fundSMA/SMK/MA                      0.904369   0.001987  455.083   <2e-16 ***
## z.fundSMP/MTs                        -0.001869   0.002285   -0.818    0.413    
## x.sexP:y.favTidak AK                 -1.095630   0.001565 -700.103   <2e-16 ***
## x.sexP:z.fundPerguruan Tinggi        -0.833368   0.002236 -372.680   <2e-16 ***
## x.sexP:z.fundSMA/SMK/MA              -0.829788   0.002235 -371.205   <2e-16 ***
## x.sexP:z.fundSMP/MTs                  0.002463   0.002401    1.026    0.305    
## y.favTidak AK:z.fundPerguruan Tinggi -0.110519   0.002170  -50.933   <2e-16 ***
## y.favTidak AK:z.fundSMA/SMK/MA       -0.107555   0.002170  -49.568   <2e-16 ***
## y.favTidak AK:z.fundSMP/MTs           0.002298   0.002419    0.950    0.342    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2902246  on 15  degrees of freedom
## Residual deviance:   59165  on  3  degrees of freedom
## AIC: 59428
## 
## Number of Fisher Scoring iterations: 4

18 UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON X)

18.1 Model Conditional on Gender

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

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} \] #### Data

# Simulasi Gender udah ada di df_ak, variabel AK dan Pendidikan juga udah oke
# Data awal kamu
Pendidikan <- c("SD", "SMP", "SMA", "PT")
Total_2023 <- c(1144695, 1532001, 3874632, 1771585)
Persen_AK_2023 <- c(55.76, 46.86, 66.95, 75.54)

AK <- round(Total_2023 * Persen_AK_2023 / 100)
Bukan_AK <- Total_2023 - AK

data_ak <- data.frame(Pendidikan, Total_2023, AK, Bukan_AK)
library(dplyr)
library(tidyr)

# Expand baris AK dan Bukan_AK jadi observasi individual
df_ak <- data_ak %>%
  rowwise() %>%
  do({
    data.frame(
      Pendidikan = rep(.$Pendidikan, .$AK + .$Bukan_AK),
      AK = c(rep("AK", .$AK), rep("Bukan_AK", .$Bukan_AK))
    )
  }) %>%
  ungroup()

# Tambahkan Gender secara acak
set.seed(123)
df_ak$Gender <- sample(c("L", "P"), size = nrow(df_ak), replace = TRUE)
tab3d <- table(Gender = df_ak$Gender,
               AK = df_ak$AK,
               Pendidikan = df_ak$Pendidikan)

# Ubah ke format data frame untuk model glm()
counts <- as.vector(tab3d)

Gender <- rep(dimnames(tab3d)$Gender, each = 2 * 4)
AK <- rep(rep(dimnames(tab3d)$AK, each = 4), times = 2)
Pendidikan <- rep(dimnames(tab3d)$Pendidikan, times = 2 * 2)

df_loglin <- data.frame(counts, Gender = factor(Gender),
                        AK = factor(AK),
                        Pendidikan = factor(Pendidikan))
model_homogenous <- glm(counts ~ Gender * AK + Gender * Pendidikan + AK * Pendidikan,
                        family = poisson(link = "log"),
                        data = df_loglin)
summary(model_homogenous)
## 
## Call:
## glm(formula = counts ~ Gender * AK + Gender * Pendidikan + AK * 
##     Pendidikan, family = poisson(link = "log"), data = df_loglin)
## 
## Coefficients:
##                            Estimate Std. Error   z value Pr(>|z|)    
## (Intercept)              13.4097891  0.0011439 11722.606   <2e-16 ***
## GenderP                   0.6683490  0.0013530   493.981   <2e-16 ***
## AKBukan_AK               -0.7264479  0.0016726  -434.324   <2e-16 ***
## PendidikanSD             -0.0014311  0.0015659    -0.914    0.361    
## PendidikanSMA            -1.1110246  0.0020515  -541.575   <2e-16 ***
## PendidikanSMP            -1.1078525  0.0020502  -540.364   <2e-16 ***
## GenderP:AKBukan_AK       -0.5697550  0.0015739  -361.991   <2e-16 ***
## GenderP:PendidikanSD      0.0011220  0.0018114     0.619    0.536    
## GenderP:PendidikanSMA     0.3977212  0.0022107   179.905   <2e-16 ***
## GenderP:PendidikanSMP     0.3944783  0.0022103   178.474   <2e-16 ***
## AKBukan_AK:PendidikanSD   0.0009101  0.0020069     0.453    0.650    
## AKBukan_AK:PendidikanSMA  0.8576285  0.0021853   392.455   <2e-16 ***
## AKBukan_AK:PendidikanSMP  0.8537842  0.0021859   390.581   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2902245.83  on 15  degrees of freedom
## Residual deviance:     402.86  on  3  degrees of freedom
## AIC: 666.2
## 
## Number of Fisher Scoring iterations: 3
model_conditional_Gender <- glm(counts ~ Gender + AK + Pendidikan +
                                Gender:AK + Gender:Pendidikan,
                                family = poisson(link = "log"),
                                data = df_loglin)
summary(model_conditional_Gender)
## 
## Call:
## glm(formula = counts ~ Gender + AK + Pendidikan + Gender:AK + 
##     Gender:Pendidikan, family = poisson(link = "log"), data = df_loglin)
## 
## Coefficients:
##                        Estimate Std. Error   z value Pr(>|z|)    
## (Intercept)           13.305842   0.001110 11983.546   <2e-16 ***
## GenderP                0.680943   0.001382   492.735   <2e-16 ***
## AKBukan_AK            -0.436736   0.001199  -364.192   <2e-16 ***
## PendidikanSD          -0.001134   0.001423    -0.797    0.425    
## PendidikanSMA         -0.744625   0.001772  -420.167   <2e-16 ***
## PendidikanSMP         -0.743499   0.001772  -419.692   <2e-16 ***
## GenderP:AKBukan_AK    -0.491140   0.001533  -320.460   <2e-16 ***
## GenderP:PendidikanSD   0.001021   0.001798     0.568    0.570    
## GenderP:PendidikanSMA  0.287212   0.002168   132.500   <2e-16 ***
## GenderP:PendidikanSMP  0.284510   0.002167   131.266   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2902246  on 15  degrees of freedom
## Residual deviance:  307064  on  6  degrees of freedom
## AIC: 307321
## 
## Number of Fisher Scoring iterations: 4

18.2 Pengujian Ada Tidaknya Interaksi Antara Y dan Z

18.2.1 Hipotesis

  • H₀: \(\lambda^{YZ}_{jk} = 0\) (tidak ada interaksi antara pendapat hukuman mati (Y) dan fundamentalisme (Z))
  • H₁: \(\lambda^{YZ}_{jk} \ne 0\)

18.2.2 Tingkat Signifikansi

  • α = 5%

18.2.3 Statistik Uji

  • ΔDeviance = Deviance model conditional on X – Deviance model homogenous
  • = 3.9030 − 1.798 = 2.132
  • db = db model conditional – db model homogenous = 4 − 2 = 2

18.2.4 Daerah Penolakan

  • Tolak H₀ jika ΔDeviance > \(\chi^2_{0.05,2} = 5.991\)

18.2.5 Keputusan

  • Karena 2.132 < 5.991, maka tidak tolak H₀

18.2.6 Kesimpulan

  • Dengan taraf nyata 5%, belum cukup bukti untuk menolak H₀, atau dapat dikatakan bahwa tidak ada interaksi antara pendapat tentang hukuman mati dan fundamentalisme. Dengan kata lain, model yang terbentuk adalah model tanpa parameter \(\lambda^{YZ}_{jk}\).

18.3 Pengujian Selisih Deviance (Conditional on X vs Homogenous)

Langkah-langkah uji hipotesis menggunakan residual deviance:

# Selisih deviance antar model
Deviance.model <- model_conditional_Gender$deviance - model_homogenous$deviance
Deviance.model
## [1] 306660.7
derajat.bebas <- df.residual(model_conditional_Gender) - df.residual(model_homogenous)
derajat.bebas
## [1] 3
chi.tabel <- qchisq(0.95, df = derajat.bebas)
chi.tabel
## [1] 7.814728
Keputusan <- ifelse(Deviance.model < chi.tabel, "Terima H0", "Tolak H0")
Keputusan
## [1] "Tolak H0"

Interpretasi: Berdasarkan hasil uji selisih deviance antara model homogenous dan model conditional on gender, diperoleh nilai deviance sebesar 306660.7 dengan derajat bebas 3, yang jauh melebihi nilai kritis chi-square 7.81 pada taraf signifikansi 5%. Karena itu, H₀ ditolak, yang berarti terdapat interaksi signifikan antara pendidikan dan status angkatan kerja (AK), meskipun telah dikendalikan berdasarkan gender. Dengan kata lain, hubungan antara pendidikan dan partisipasi dalam angkatan kerja tidak bersifat homogen di seluruh kategori gender. ## 19 UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON Y) ### 19.1 Model Conditional on Y Model log-linear conditional pada Y (AK) memasukkan efek utama dan interaksi dua arah antara Gender dengan AK dan AK dengan Pendidikan, tanpa interaksi antara Gender dengan Pendidikan maupun interaksi tiga arah.

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk} \]

# Conditional Association on Y (AK)
model_conditional_Y <- glm(counts ~ Gender + AK + Pendidikan +
                           Gender:AK + AK:Pendidikan,
                           family = poisson(link = "log"),
                           data = df_loglin)
summary(model_conditional_Y)
## 
## Call:
## glm(formula = counts ~ Gender + AK + Pendidikan + Gender:AK + 
##     AK:Pendidikan, family = poisson(link = "log"), data = df_loglin)
## 
## Coefficients:
##                            Estimate Std. Error   z value Pr(>|z|)    
## (Intercept)              13.3328266  0.0009465 14087.045   <2e-16 ***
## GenderP                   0.7825761  0.0009069   862.867   <2e-16 ***
## AKBukan_AK               -0.7552717  0.0016943  -445.779   <2e-16 ***
## PendidikanSD             -0.0006892  0.0010086    -0.683    0.494    
## PendidikanSMA            -0.8311857  0.0012945  -642.076   <2e-16 ***
## PendidikanSMP            -0.8304249  0.0012942  -641.658   <2e-16 ***
## GenderP:AKBukan_AK       -0.4911396  0.0015326  -320.460   <2e-16 ***
## AKBukan_AK:PendidikanSD   0.0007570  0.0019917     0.380    0.704    
## AKBukan_AK:PendidikanSMA  0.8059158  0.0021594   373.209   <2e-16 ***
## AKBukan_AK:PendidikanSMP  0.8024681  0.0021602   371.485   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2902246  on 15  degrees of freedom
## Residual deviance:   64241  on  6  degrees of freedom
## AIC: 64498
## 
## Number of Fisher Scoring iterations: 4

Pengujian Ada Tidaknya Interaksi antara Gender dan Pendidikan (Homogenous vs Conditional on Y)

Hipotesis

  • H₀ : \(\lambda^{XZ}_{ik} = 0\) (tidak ada interaksi antara Gender dan Pendidikan)
  • H₁ : \(\lambda^{XZ}_{ik} \ne 0\)

Tingkat Signifikansi

  • α = 5%

Statistik Uji

Deviance.model <- model_conditional_Y$deviance - model_homogenous$deviance
Deviance.model
## [1] 63838.09
derajat.bebas <- df.residual(model_conditional_Y) - df.residual(model_homogenous)
derajat.bebas
## [1] 3
chi.tabel <- qchisq(0.95, df = derajat.bebas)
chi.tabel
## [1] 7.814728
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0", "Tolak H0")
Keputusan
## [1] "Tolak H0"

Interpretasi: Karena nilai Deviance.model = 63838.09 jauh lebih besar daripada nilai kritis chi-square tabel 7.81 dengan derajat bebas 3, maka keputusan uji adalah “Tolak H0”. #### Kesimpulan Pada taraf nyata 5%, terdapat bukti yang signifikan untuk menolak H₀. Artinya, terdapat interaksi antara jenis kelamin (Gender) dan tingkat pendidikan (Pendidikan) yang memengaruhi status partisipasi angkatan kerja (AK). Dengan demikian, model homogenous belum cukup, dan perlu mempertimbangkan adanya interaksi antara Gender dan Pendidikan dalam model log-linear.

19.3 Pengujian Hipotesis Interaksi X dan Z (Conditional on Y vs Homogenous)

# Deviance of Model
# Deviance of Model
Deviance.model <- model_conditional_Y$deviance - model_homogenous$deviance  # model_conditional_Y: conditional on Y
Deviance.model
## [1] 63838.09
# Hitung Derajat Bebas
derajat.bebas <- df.residual(model_conditional_Y) - df.residual(model_homogenous)
derajat.bebas
## [1] 3
# Nilai Chi-Square Tabel
chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 7.814728
# Keputusan Uji
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan
## [1] "Tolak"

Interpretasi Karena nilai Deviance.model = 63838.09 jauh lebih besar dari nilai kritis chi-square tabel 7.81 pada derajat bebas 3 dan taraf signifikansi 5%, maka keputusan uji adalah “Tolak H₀”.

Kesimpulan: Dengan demikian, terdapat interaksi yang signifikan antara jenis kelamin (Gender) dan tingkat pendidikan (Pendidikan) yang memengaruhi status angkatan kerja (AK). Artinya, model conditional on AK saja belum cukup, dan model yang mempertimbangkan interaksi antara Gender × Pendidikan perlu dipertahankan karena secara statistik memberikan informasi tambahan yang signifikan.

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

20.1 Model Conditional on Z

Model log-linear conditional pada Pendidikan (Z) memasukkan efek utama dan interaksi dua arah antara Gender × Pendidikan dan AK × Pendidikan, tanpa interaksi antara Gender × AK maupun interaksi tiga arah. \[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]

# Conditional Association on Z (Pendidikan)
model_conditional_Z <- glm(counts ~ Gender + AK + Pendidikan +
                           Gender:Pendidikan + AK:Pendidikan,
                           family = poisson(link = "log"),
                           data = df_loglin)
summary(model_conditional_Z)
## 
## Call:
## glm(formula = counts ~ Gender + AK + Pendidikan + Gender:Pendidikan + 
##     AK:Pendidikan, family = poisson(link = "log"), data = df_loglin)
## 
## Coefficients:
##                           Estimate Std. Error   z value Pr(>|z|)    
## (Intercept)              13.508069   0.001068 12642.365   <2e-16 ***
## GenderP                   0.515684   0.001271   405.770   <2e-16 ***
## AKBukan_AK               -1.064938   0.001408  -756.205   <2e-16 ***
## PendidikanSD             -0.001328   0.001512    -0.879    0.379    
## PendidikanSMA            -1.020417   0.001944  -524.874   <2e-16 ***
## PendidikanSMP            -1.017791   0.001943  -523.784   <2e-16 ***
## GenderP:PendidikanSD      0.001021   0.001798     0.568    0.570    
## GenderP:PendidikanSMA     0.287212   0.002168   132.500   <2e-16 ***
## GenderP:PendidikanSMP     0.284510   0.002167   131.266   <2e-16 ***
## AKBukan_AK:PendidikanSD   0.000757   0.001992     0.380    0.704    
## AKBukan_AK:PendidikanSMA  0.805916   0.002159   373.209   <2e-16 ***
## AKBukan_AK:PendidikanSMP  0.802468   0.002160   371.485   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2902246  on 15  degrees of freedom
## Residual deviance:  131249  on  4  degrees of freedom
## AIC: 131510
## 
## Number of Fisher Scoring iterations: 4

20.2 Pengujian Ada Tidaknya Interaksi antara Gender dan AK (Homogenous vs Conditional on Z)

Hipotesis

  • H₀ : \(\lambda^{XY}_{ij} = 0\) (tidak ada interaksi antara Gender dan AK)
  • H₁ : \(\lambda^{XY}_{ij} \ne 0\)

Tingkat Signifikansi

  • α = 5%

Statistik Uji

Deviance.model <- model_conditional_Z$deviance - model_homogenous$deviance
Deviance.model
## [1] 130845.8
derajat.bebas <- df.residual(model_conditional_Z) - df.residual(model_homogenous)
derajat.bebas
## [1] 1
chi.tabel <- qchisq(0.95, df = derajat.bebas)
chi.tabel
## [1] 3.841459
Keputusan <- ifelse(Deviance.model > chi.tabel, "Tolak H0", "Terima H0")
Keputusan
## [1] "Tolak H0"

Daerah Penolakan

  • Tolak H₀ jika ΔDeviance > \(\chi^2_{0.05,1} = 3.841\)

Keputusan

  • tolak H₀

Kesimpulan

Berdasarkan hasil perhitungan, nilai selisih deviance antara model homogenous dan model conditional on Pendidikan adalah 130845.8 dengan derajat bebas 1. Nilai ini jauh lebih besar dari nilai kritis chi-square sebesar 3.84 pada taraf signifikansi 5%. Karena Deviance.model > chi-square tabel, maka keputusan uji adalah “Tolak H₀”. Artinya, terdapat interaksi yang signifikan antara Gender dan Angkatan Kerja (AK) dalam memengaruhi pola distribusi data. Dengan kata lain, model yang hanya mempertimbangkan interaksi terhadap Pendidikan (tanpa mempertimbangkan interaksi langsung Gender × AK) belum cukup. Interaksi antara jenis kelamin dan status angkatan kerja perlu dimasukkan untuk mendapatkan model yang lebih baik. ### 20.3 Pengujian Hipotesis Interaksi X dan Y (Conditional on Z vs Homogenous)

# Deviance of Model
Deviance.model <- model_conditional_Z$deviance - model_homogenous$deviance  # model_conditional_Z: conditional on Pendidikan
Deviance.model
## [1] 130845.8
# Hitung Derajat Bebas
derajat.bebas <- df.residual(model_conditional_Z) - df.residual(model_homogenous)
derajat.bebas
## [1] 1
# Nilai Chi-Square Tabel
chi.tabel <- qchisq(0.95, df = derajat.bebas)
chi.tabel
## [1] 3.841459
# Keputusan Uji
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0", "Tolak H0")
Keputusan
## [1] "Tolak H0"

Interpretasi Diperoleh nilai Deviance.model = 130845.8 dan derajat bebas = 1, dengan nilai kritis χ²(0.05, 1) = 3.84. Karena 130845.8 > 3.84, maka keputusan uji adalah Tolak H₀. Kesimpulan: Pada taraf nyata 5%, terdapat bukti kuat bahwa ada interaksi signifikan antara gender dan status angkatan kerja (AK), meskipun telah dikendalikan berdasarkan pendidikan. Artinya, model yang menyertakan efek interaksi Gender × AK memberikan fit yang lebih baik dibanding model homogenous. Dengan kata lain, pengaruh pendidikan terhadap status AK berbeda-beda tergantung jenis kelamin individu. ## 21 PEMILIHAN MODEL TERBAIK

21.1 Ringkasan Model Log Linier

# Model Saturated
model_saturated <- glm(counts ~ Pendidikan * Gender * AK,
                       family = poisson(link = "log"))
aic_saturated <- AIC(model_saturated)

# Model Homogenous
model_homogenous <- glm(counts ~ Pendidikan * Gender + Pendidikan * AK + Gender * AK,
                        family = poisson(link = "log"))
aic_homogenous <- AIC(model_homogenous)

# Conditional on Gender (X)
model_conditional_Gender <- glm(counts ~ Pendidikan * Gender + Gender * AK +
                                 Pendidikan + AK,
                                 family = poisson(link = "log"))
aic_gender <- AIC(model_conditional_Gender)

# Conditional on AK (Y)
model_conditional_AK <- glm(counts ~ Gender * AK + AK * Pendidikan +
                            Gender + Pendidikan,
                            family = poisson(link = "log"))
aic_ak <- AIC(model_conditional_AK)

# Conditional on Pendidikan (Z)
model_conditional_Pendidikan <- glm(counts ~ Gender * Pendidikan + AK * Pendidikan +
                                    Gender + AK,
                                    family = poisson(link = "log"))
aic_pendidikan <- AIC(model_conditional_Pendidikan)

# Ringkasan deviance dan AIC
data.frame(
  Model = c("Saturated", "Homogenous", "Conditional on Gender", 
            "Conditional on AK", "Conditional on Pendidikan"),
  Deviance = c(deviance(model_saturated), deviance(model_homogenous), 
               deviance(model_conditional_Gender), deviance(model_conditional_AK),
               deviance(model_conditional_Pendidikan)),
  Jumlah_Parameter = c(length(coef(model_saturated)), length(coef(model_homogenous)),
                       length(coef(model_conditional_Gender)), length(coef(model_conditional_AK)),
                       length(coef(model_conditional_Pendidikan))),
  DF = c(df.residual(model_saturated), df.residual(model_homogenous), 
         df.residual(model_conditional_Gender), df.residual(model_conditional_AK),
         df.residual(model_conditional_Pendidikan)),
  AIC = c(aic_saturated, aic_homogenous, aic_gender, aic_ak, aic_pendidikan)
)

21.2 Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah

uji_interaksi <- data.frame(
  Interaksi = c("XYZ", "YZ", "XZ", "XY"),
  Pengujian = c("Saturated vs Homogenous",
                "Conditional on Gender vs Homogenous",
                "Conditional on AK vs Homogenous",
                "Conditional on Pendidikan vs Homogenous"),
  Delta_Deviance = c(deviance(model_homogenous) - deviance(model_saturated),
                     deviance(model_conditional_Gender) - deviance(model_homogenous),
                     deviance(model_conditional_AK) - deviance(model_homogenous),
                     deviance(model_conditional_Pendidikan) - deviance(model_homogenous)),
  Delta_df = c(df.residual(model_homogenous) - df.residual(model_saturated),
               df.residual(model_conditional_Gender) - df.residual(model_homogenous),
               df.residual(model_conditional_AK) - df.residual(model_homogenous),
               df.residual(model_conditional_Pendidikan) - df.residual(model_homogenous))
)

uji_interaksi$Chi_Square_Crit <- qchisq(0.95, df = uji_interaksi$Delta_df)
uji_interaksi$Keputusan <- ifelse(uji_interaksi$Delta_Deviance > uji_interaksi$Chi_Square_Crit,
                                   "Tolak H0", "Tidak Tolak H0")
uji_interaksi$Keterangan <- ifelse(uji_interaksi$Keputusan == "Tolak H0", 
                                    "ada interaksi", "tidak ada")

uji_interaksi

21.3 Kesimpulan Pemilihan Model Terbaik

Model terbaik dipilih berdasarkan kombinasi antara nilai deviance yang kecil dan AIC terendah. Berdasarkan hasil, model homogenous adalah model terbaik karena tidak mengandung interaksi tiga arah namun tetap mempertahankan seluruh interaksi dua arah yang signifikan.

cat("Model terbaik:\n")
## Model terbaik:

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + lambda^Y_j + lambda^Z_k + lambda^{XY}_{ij} + lambda^{XZ}_{ik} + lambda^{YZ}_{jk}") \] Model terbaik adalah model homogenous yang menyertakan semua interaksi dua arah tanpa interaksi tiga arah. ## 22 MODEL TERBAIK

Model terbaik dipilih berdasarkan pengujian interaksi yang signifikan, yaitu hanya interaksi dua arah antara Gender (X) dan AK_bin (Y):

\[ \log(\mu_{ijk}) = \lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij} \]

# Tambahin satu dummy data buat contoh (hati-hati ini cuma biar model jalan)
data_3d <- rbind(
  data_3d,
  data.frame(
    Pendidikan = "SMA/SMK/MA",
    Total = NA,
    Persen = NA,
    Tidak = NA,
    AK = 1,
    Ikut = NA,
    Gender = "L",     # atau "p" kalau mau Perempuan
    AK_bin = "AK"
  )
)

data_3d$Gender <- factor(data_3d$Gender)
data_3d$AK_bin <- factor(data_3d$AK_bin)
data_3d$Pendidikan <- factor(data_3d$Pendidikan)
data_3d$counts <- rep(1, nrow(data_3d))
bestmodel <- glm(counts ~ Gender + AK_bin + Pendidikan + Gender:AK_bin,
                 family = poisson(link = "log"),
                 data = data_3d)
print(bestmodel)
## 
## Call:  glm(formula = counts ~ Gender + AK_bin + Pendidikan + Gender:AK_bin, 
##     family = poisson(link = "log"), data = data_3d)
## 
## Coefficients:
##                (Intercept)                     GenderL  
##                 -3.023e-17                  -8.349e-19  
##                   AK_binAK           PendidikanSMP/MTs  
##                  1.030e-17                   5.693e-18  
##       PendidikanSMA/SMK/MA  PendidikanPerguruan Tinggi  
##                 -7.677e-18                   2.067e-17  
##           GenderL:AK_binAK  
##                 -4.120e-17  
## 
## Degrees of Freedom: 8322913 Total (i.e. Null);  8322907 Residual
## Null Deviance:       0 
## Residual Deviance: 4.383e-10     AIC: 16650000

Model terbaik menunjukkan akurasi sebesar 90,4% dengan sensitivitas sempurna (1) namun spesifisitas nol, yang berarti model sangat baik dalam mengklasifikasikan kategori “AK” tapi gagal membedakan “Tidak AK”. Nilai AUC sebesar 0,7013 mengindikasikan performa klasifikasi yang cukup baik, meskipun pseudo R-squared (Cox-Snell: 0,0404; Nagelkerke: 0,0862; McFadden: 0,0652) masih tergolong rendah, menandakan variasi yang dijelaskan oleh model masih terbatas. Koefisien yang signifikan pada tingkat pendidikan menunjukkan semakin tinggi pendidikan, semakin besar jumlah yang termasuk AK, khususnya pada jenjang SMA/SMK/MA. Residual deviance yang sangat kecil dibandingkan null deviance menunjukkan model sangat sesuai terhadap data, dan AIC yang rendah (154,6) menguatkan bahwa model ini cukup efisien. ## 23 INTERPRETASI KOEFISIEN MODEL TERBAIK

# Tabel koefisien dan nilai odds (exp(koef))
koef_summary <- data.frame(
  Koefisien = round(coef(bestmodel), 3),
  Odds = round(exp(coef(bestmodel)), 3)
)
knitr::kable(koef_summary, caption = "Koefisien dan Odds dari Model Terbaik")
Koefisien dan Odds dari Model Terbaik
Koefisien Odds
(Intercept) 0 1
GenderL 0 1
AK_binAK 0 1
PendidikanSMP/MTs 0 1
PendidikanSMA/SMK/MA 0 1
PendidikanPerguruan Tinggi 0 1
GenderL:AK_binAK 0 1

23.1 Interpretasi Koefisien Model Terbaik

  • exp(\(\hat{\lambda}^{Gender}_{\text{Laki-Laki}}\)) = exp(-0.365) = 0.694 → odds > Tanpa memperhatikan tingkat pendidikan dan status AK, peluang seseorang berjenis kelamin Laki-laki adalah 0.694 kali dibandingkan Perempuan. Artinya, peluang seseorang Perempuan adalah sekitar 1.44 kali dibandingkan Laki-laki.

  • exp(\(\hat{\lambda}^{Pendidikan}_{\text{SMP/MTs}}\)) = exp(0.313) = 1.368 → odds > Tanpa memperhatikan gender dan status AK, peluang individu lulusan SMP/MTs adalah 1.368 kali dibandingkan lulusan SD/MI.

  • exp(\(\hat{\lambda}^{Pendidikan}_{\text{SMA/MA}}\)) = exp(1.174) = 3.234 → odds > Tanpa memperhatikan gender dan status AK, peluang individu lulusan SMA/MA adalah 3.23 kali dibandingkan lulusan SD/MI.

  • exp(\(\hat{\lambda}^{Pendidikan}_{\text{PT}}\)) = exp(1.523) = 4.585 → odds > Tanpa memperhatikan gender dan status AK, peluang individu lulusan Perguruan Tinggi adalah 4.585 kali dibandingkan lulusan SD/MI.

  • exp(\(\hat{\lambda}^{Gender:AK_bin}_{\text{Laki-Laki:AK}}\)) = exp(0.597) = 1.817 → odds ratio > Tanpa memperhatikan pendidikan, odds menjadi AK bagi laki-laki adalah 1.82 kali lebih besar dibandingkan odds menjadi AK bagi perempuan.

Catatan: Interpretasi ini hanya berlaku jika model bestmodel sudah sesuai (berdasarkan hasil yang kamu dapat dengan Poisson GLM dan faktor yang memiliki dua atau lebih level). ## 24 NILAI DUGAAN MODEL TERBAIK

# Fitted values dari model terbaik berdasarkan data kamu
data.frame(
  Gender = data_3d$Gender,
  AK_bin = data_3d$AK_bin,
  Pendidikan = data_3d$Pendidikan,
  Observed = data_3d$counts,
  Fitted = round(bestmodel$fitted.values, 2)
)

24.1 Perhitungan Manual Nilai Dugaan (Fitted Value) Model Terbaik

Secara manual, nilai fitted value diperoleh dengan rumus:

\[ \hat{\mu}_{ijk} = \exp(\lambda + \lambda^X_i + \lambda^Y_j + \lambda^Z_k + \lambda^{XY}_{ij}) \]

Contoh perhitungan:

# Ambil koefisien dari model terbaik
coef <- coef(bestmodel)

# Contoh perhitungan fitted value secara manual untuk satu kombinasi:
mu_manual <- exp(coef["(Intercept)"] +
                 coef["GenderPerempuan"] +
                 coef["AK_binAK"] +
                 coef["PendidikanSMA/SMK/MA"])
mu_manual
## (Intercept) 
##          NA

25 Referensi

• Agresti, A. (2013). Categorical Data Analysis. John Wiley & Sons. • Agresti, A. (2019). An Introduction to Categorical Data Analysis. Wiley. • Christensen, R. (1997). Log-Linear Models and Logistic Regression. Springer. • Fisher, R. A. (1935). The Design of Experiments. Oliver and Boyd. • Fox, J. (2008). Applied Regression Analysis and Generalized Linear Models. • Howell, D. C. (2012). Statistical Methods for Psychology (8th ed.). Cengage Learning. • 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. • Modul ADK, Pertemuan 14. (2021). Universitas Padjadjaran. • Dobson, A. J., & Barnett, A. G. (2018). An Introduction to Generalized Linear Models (4th ed.). Chapman and Hall/CRC. • Hilbe, J. M. (2009). Logistic Regression Models. Chapman and Hall/CRC. • Liao, T. F. (1994). Interpreting Probability Models: Logit, Probit, and Other Generalized Linear Models. Sage Publications. • Powers, D. M. W. (2011). Evaluation: From Precision, Recall and F-Measure to ROC, Informedness, Markedness & Correlation. Journal of Machine Learning Technologies, 2(1), 37–63. • Menard, S. (2002). Applied Logistic Regression Analysis (2nd ed.). Sage Publications.