Pendahuluan

Dalam era informasi seperti saat ini, data kategori menjadi bagian penting dalam berbagai bidang, seperti kesehatan, ekonomi, pendidikan, dan ilmu sosial. Data jenis ini muncul ketika informasi diklasifikasikan ke dalam kategori atau kelompok, bukan dalam bentuk angka kontinu. Contoh sederhananya adalah jenis kelamin (laki-laki/perempuan), status merokok (ya/tidak), atau tingkat kepuasan (rendah/sedang/tinggi).

Analisis data kategori berperan penting dalam memahami hubungan antar kategori, menguji asumsi independensi, serta mengukur kekuatan asosiasi di antara variabel. Tidak seperti data numerik yang dapat langsung dianalisis dengan teknik statistik parametrik, data kategori memerlukan pendekatan khusus seperti uji proporsi, analisis tabel kontingensi, dan model regresi khusus seperti Generalized Linear Models (GLM).

Selain itu, seiring kemajuan teknologi dan era big data, teknik analisis data kategori juga semakin terintegrasi dalam aplikasi machine learning seperti decision tree, random forest, dan analisis korespondensi. Artinya, pemahaman terhadap data kategori menjadi semakin penting dalam praktik analisis data modern, baik untuk eksplorasi data maupun pembangunan model prediktif.

Tujuan Analisis Data Kategori

  1. Mengungkap pola tersembunyi dalam data kategori
    Dengan klasifikasi dan eksplorasi kategori, kita bisa menemukan tren atau hubungan tertentu. Misalnya, analisis preferensi belanja menunjukkan bahwa perempuan lebih cenderung memilih belanja online pada hari kerja.

  2. Menilai hubungan antarvariabel kategori
    Apakah perokok lebih cenderung mengalami hipertensi dibandingkan non-perokok? Dengan analisis tabel kontingensi dan uji asosiasi, kita dapat menilai kekuatan hubungan tersebut.

  3. Menyediakan dasar bagi pengambilan keputusan
    Data kategori sering digunakan dalam survei publik. Hasil dari variabel seperti “setuju/tidak setuju” bisa menjadi dasar kebijakan pemerintah dalam merancang program sosial.

  4. Membangun model prediktif untuk kejadian tertentu
    Seperti regresi logistik untuk memprediksi kemungkinan seorang pelanggan akan melakukan pembelian atau tidak berdasarkan kategori-kategori seperti jenis kelamin, penghasilan, atau domisili.

Definisi dan Ruang Lingkup

Analisis data kategori berfokus pada variabel yang nilainya berupa kategori, baik nominal maupun ordinal.

  • Nominal: Kategori tidak memiliki urutan.
    Contoh: warna favorit (merah, biru, kuning), jenis kelamin (L/P).

  • Ordinal: Kategori memiliki urutan logis.
    Contoh: tingkat pendidikan (SD < SMP < SMA < S1), atau kepuasan pelanggan (rendah < sedang < tinggi).

Jenis Data Kategori

  • Biner
    Memiliki dua nilai kategori.
    Contoh: ya/tidak, hadir/tidak hadir, positif/negatif. Biasa dianalisis dengan regresi logistik biner.

  • Multikategori
    Memiliki tiga atau lebih nilai kategori.
    Contoh: status ekonomi (rendah/sedang/tinggi), tingkat pendidikan.
    Untuk ini, digunakan regresi logistik multinomial atau model log-linear.

Perbedaan dengan Data Kuantitatif

Aspek Data Kuantitatif Data Kategorik
Nilai Angka kontinu/diskret Kategori (label)
Metode Analisis Parametrik: t-test, ANOVA, regresi Non-parametrik: chi-square, GLM
Visualisasi Histogram, boxplot, scatterplot Bar chart, mosaic plot
Contoh Usia, tinggi badan, pendapatan Jenis kelamin, status pekerjaan

Data kuantitatif bisa dirata-ratakan, tapi data kategori hanya bisa dihitung frekuensinya.

Manfaat Analisis Data Kategori dalam Berbagai Bidang

Ilmu Sosial & Psikologi

  • Menilai hasil survei: seberapa banyak responden setuju dengan suatu isu.
  • Melakukan segmentasi sosial: seperti berdasarkan tingkat kepercayaan terhadap pemerintah.

Kesehatan & Medis

  • Mengelompokkan pasien berdasarkan gejala (ringan/sedang/berat).
  • Menguji hubungan antara perilaku (merokok) dengan penyakit (kanker paru).

Bisnis & Pemasaran

  • Menganalisis loyalitas pelanggan dari kategori umur dan gender.
  • Mengembangkan strategi pemasaran berdasarkan tipe konsumen (baru/lama/potensial).

Pendidikan

  • Menganalisis metode pembelajaran berdasarkan prestasi (rendah/sedang/tinggi).
  • Menilai efektivitas program bimbingan berdasarkan latar belakang siswa.

Pemerintahan

  • Menyusun program bantuan sosial berdasar kategori pendapatan warga.
  • Membuat kebijakan publik dari hasil polling masyarakat (setuju/tidak).

Keamanan

  • Analisis kejahatan berdasarkan wilayah (urban/suburban) dan waktu kejadian.
  • Klasifikasi pelaku kejahatan berdasarkan status sosial atau jenis kelamin.

Metode Umum dalam Analisis Data Kategori

  1. Tabel Kontingensi & Uji Chi-Square
    Dasar dari semua analisis kategori. Digunakan untuk menguji apakah dua variabel kategori saling bebas atau tidak.
    Contoh: apakah jenis kelamin berhubungan dengan pilihan jurusan?

  2. Regresi Logistik
    Cocok untuk prediksi kejadian biner seperti “lulus/tidak lulus”, “membeli/tidak membeli”. Menghasilkan interpretasi dalam bentuk peluang dan odds ratio.

  3. Correspondence Analysis (CA)
    Metode eksploratif untuk melihat hubungan visual antara baris dan kolom tabel kontingensi. Cocok untuk banyak kategori.

  4. Decision Tree & Random Forest
    Metode machine learning berbasis pemisahan data secara berulang. Sangat efektif untuk klasifikasi berdasarkan variabel kategori dan kombinasi kategorik-kuantitatif.

Sampling Design dan Ukuran Asosiasi

Bagian ini membahas bagaimana desain pengumpulan data memengaruhi jenis analisis yang bisa dilakukan serta bagaimana menilai kekuatan hubungan antar dua variabel kategori.

Desain Sampling

Desain sampling adalah strategi dalam pengumpulan data yang akan memengaruhi validitas inferensi statistik. Dua desain yang paling umum:

Desain Prospektif (Forward-looking)

  • Data dikumpulkan dengan mengikuti subjek ke depan waktu.
  • Contoh: studi kohort atau eksperimen klinis.
  • Contoh Kasus: Peneliti mengelompokkan subjek berdasarkan status merokok, lalu mengikuti selama 5 tahun untuk mencatat apakah mereka terkena penyakit jantung.

Karakteristik:

  • Dapat menghitung probabilitas kejadian secara langsung.
  • Ukuran asosiasi: RD dan RR.

Desain Retrospektif (Backward-looking)

  • Dimulai dari outcome, kemudian menelusuri faktor risikonya.
  • Contoh: studi kasus-kontrol.
  • Contoh Kasus: Peneliti mengambil pasien kanker paru dan mencocokkan dengan kontrol sehat, lalu menanyakan riwayat merokok.

Karakteristik: - Tidak bisa menghitung probabilitas populasi. - Ukuran yang digunakan: OR.

Ukuran Asosiasi

Ukuran asosiasi digunakan untuk mengukur kekuatan dan arah hubungan antara dua variabel dalam tabel 2x2. Tiga ukuran utama:

Risk Difference (RD)

Risk Difference (RD) atau selisih risiko mengukur perbedaan probabilitas kejadian antara dua kelompok. RD dihitung sebagai:

\[ RD = P(kejadian | grup\ 1) - P(kejadian | grup\ 2) \]

Interpretasi:

  • RD > 0: Risiko kejadian lebih tinggi di grup 1.
  • RD < 0: Risiko kejadian lebih tinggi di grup 2.
  • RD = 0: Tidak ada perbedaan risiko antar kelompok.
RD <- function(n11, n12, n21, n22) {
  (n11 / (n11 + n12)) - (n21 / (n21 + n22))
}
RD(50, 50, 30, 70)
## [1] 0.2

Relative Risk (RR)

Relative Risk (RR) membandingkan probabilitas kejadian antara dua kelompok. Ukuran ini sangat umum dalam studi kohort dan eksperimental.

\[ RR = \frac{P(kejadian | grup\ 1)}{P(kejadian | grup\ 2)} \]

Interpretasi: - RR > 1: Risiko di grup 1 lebih tinggi. - RR < 1: Risiko di grup 2 lebih tinggi. - RR = 1: Risiko kejadian sama di kedua grup.

RR <- function(n11, n12, n21, n22) {
  (n11 / (n11 + n12)) / (n21 / (n21 + n22))
}
RR(50, 50, 30, 70)
## [1] 1.666667

Odds Ratio (OR)

Odds Ratio (OR) mengukur perbandingan odds (peluang kejadian dibandingkan dengan tidak kejadian) antara dua kelompok. Umumnya digunakan dalam studi kasus-kontrol.

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

Interpretasi:

  • OR > 1: Odds kejadian lebih tinggi di grup 1. - OR < 1: Odds kejadian lebih tinggi di grup 2. - OR = 1: Odds kejadian sama.

Contoh Kasus:

Sebuah studi kasus-kontrol ingin mengetahui apakah ada hubungan antara kebiasaan merokok dan kejadian penyakit paru-paru.

Berikut ini adalah data hasil penelitian terhadap 200 orang:

Penyakit Paru (+) Tidak Sakit (-) Total
Perokok 50 50 100
Non-Perokok 30 70 100 
OR <- function(n11, n12, n21, n22) {
  (n11 * n22) / (n12 * n21)
}
OR(50, 50, 30, 70)
## [1] 2.333333

Kesimpulan:

Karena OR > 1, odds terjadinya penyakit paru-paru lebih tinggi pada kelompok perokok dibandingkan non-perokok. Artinya, merokok memiliki asosiasi positif terhadap risiko penyakit paru-paru, dan dapat dianggap sebagai faktor risiko berdasarkan data ini.

Inferensi dalam Analisis Data Kategori

Estimasi Proporsi

  • Estimasi titik dan interval

Uji Proporsi, Asosiasi, dan Independensi

Contoh kode uji proporsi atau chi-square

data <- matrix(c(30, 10, 15, 45), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Terpapar" = c("Ya", "Tidak"), "Kejadian" = c("Ya", "Tidak"))
chisq.test(data)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data
## X-squared = 22.264, df = 1, p-value = 2.376e-06

Hasil dan Interpretasi

Kejadian = Ya Kejadian = Tidak
Terpapar = Ya 30 10
Terpapar = Tidak 15 45
  • Statistik chi-square: 22.264
  • df: 1
  • p-value: 2.376e-06

Kesimpulan:

Karena p-value < 0.05, maka tolak H0. Artinya, ada hubungan signifikan antara paparan dan kejadian. Orang yang terpapar cenderung lebih sering mengalami kejadian dibanding yang tidak terpapar.

Distribusi Probabilitas dalam Data Kategori

Variabel acak kategorik adalah variabel yang nilainya terbatas pada sejumlah kategori yang bersifat diskrit. Untuk memahami bagaimana peluang tiap kategori muncul, digunakan distribusi probabilitas yang sesuai dengan bentuk data kategorik tersebut. Pada bagian ini, akan dibahas empat distribusi utama yang sering digunakan dalam analisis data kategori, yaitu distribusi Bernoulli, Binomial, Multinomial, dan Poisson.

Distribusi Bernoulli

Distribusi Bernoulli digunakan untuk menggambarkan percobaan yang hanya memiliki dua hasil yang mungkin: “sukses” (1) dan “gagal” (0). Probabilitas terjadinya sukses dilambangkan dengan \(p\), sedangkan kegagalan sebesar \(1 - p\).

Fungsi probabilitas:

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

Contoh kasus:

Seorang siswa mengerjakan satu soal pilihan ganda dengan dua opsi jawaban. Distribusi Bernoulli dapat digunakan untuk memodelkan peluang siswa menjawab dengan benar.

library(knitr)
library(kableExtra)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
## 
##     group_rows
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
set.seed(123)
bernoulli_sample <- rbinom(n = 10, size = 1, prob = 0.7)
bernoulli_sample
##  [1] 1 0 1 0 0 1 1 0 1 1

Kesimpulan:

Dari 10 kali simulasi percobaan:

Nilai 1 muncul sebanyak 6 kali → menandakan siswa menjawab benar.

Nilai 0 muncul sebanyak 4 kali → menandakan siswa menjawab salah.

Karena probabilitas menjawab benar adalah 0.7, hasil simulasi ini konsisten dengan ekspektasi, yaitu bahwa siswa lebih sering menjawab dengan benar dibandingkan salah. Meskipun hasil simulasi tidak selalu menghasilkan persentase yang sama persis dengan probabilitas teoretis, secara umum pola distribusinya mencerminkan bahwa siswa cenderung memiliki peluang tinggi untuk menjawab dengan benar.

Distribusi Binomial

Distribusi Binomial adalah generalisasi dari distribusi Bernoulli yang digunakan untuk \(n\) kali percobaan independen yang masing-masing memiliki dua kemungkinan hasil (sukses atau gagal). Jika setiap percobaan memiliki probabilitas sukses \(p\), maka distribusi Binomial memiliki fungsi probabilitas:

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

Keterangan Notasi:

  • \(X\): Jumlah keberhasilan dalam \(n\) percobaan
  • \(n\): Jumlah percobaan
  • \(k\): Jumlah keberhasilan yang diamati
  • \(p\): Probabilitas keberhasilan dalam satu percobaan
  • \(\binom{n}{k}\): Kombinasi “n pilih k”, dihitung sebagai:
    \[ \binom{n}{k} = \frac{n!}{k!(n-k)!} \]

Contoh Variabel Acak Binomial - Jumlah keberhasilan dalam 10 kali lemparan koin. - Jumlah pasien yang sembuh setelah diberikan obat tertentu dalam suatu studi klinis.

Contoh Kasus

Seorang peneliti ingin mengetahui berapa banyak dari 8 siswa yang akan lulus ujian, jika probabilitas lulus tiap siswa adalah 0.6. Simulasikan 10 kali percobaan dengan distribusi binomial.

set.seed(123)
binomial_sample <- rbinom(n = 10, size = 8, prob = 0.6)
binomial_sample
##  [1] 6 4 5 3 3 7 5 3 5 5

Kesimpulan:

Dari 10 kali simulasi:

Nilai lulus paling sering muncul adalah 5 siswa

Jumlah siswa yang lulus berkisar antara 3 hingga 7 orang

Hal ini menunjukkan bahwa simulasi dengan distribusi binomial mencerminkan kemungkinan nyata bahwa rata-rata 5 dari 8 siswa akan lulus jika peluang kelulusan masing-masing adalah 0.6. Karena proses bersifat acak, hasilnya bisa bervariasi, tapi akan cenderung mendekati nilai harapan dalam banyak percobaan.

Distribusi Multinomial

Distribusi Multinomial merupakan perluasan dari distribusi Binomial yang digunakan ketika suatu percobaan memiliki lebih dari dua kemungkinan hasil. Jika sebuah eksperimen dilakukan sebanyak \(n\) kali, dan setiap percobaan dapat menghasilkan salah satu dari \(k\) kategori dengan probabilitas masing-masing \(p_1, p_2, \dots, p_k\), maka distribusi ini menjadi model yang tepat untuk menggambarkan peluang kombinasi hasil dari keseluruhan percobaan tersebut.

Fungsi probabilitas:

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

  • \(x_i\): jumlah kejadian kategori ke-\(i\)
  • \(p_i\): probabilitas kategori ke-\(i\)
  • \(\sum x_i = n\), dan \(\sum p_i = 1\)

Contoh kasus :

Sebuah kedai kopi mencatat pilihan minuman dari 15 pelanggan. Terdapat tiga jenis minuman yang ditawarkan: kopi, teh, dan cokelat. Probabilitas masing-masing minuman diperkirakan sebagai berikut:

Kopi: 40%

Teh: 40%

Cokelat: 20%

Distribusi multinomial digunakan untuk memodelkan kemungkinan jumlah pelanggan yang memilih tiap jenis minuman.

set.seed(123)
multinomial_sample <- rmultinom(n = 1, size = 15, prob = c(0.4, 0.4, 0.2))
multinomial_sample
##      [,1]
## [1,]    5
## [2,]    5
## [3,]    5

Kesimpulan

Hasil simulasi menunjukkan bahwa dari 15 pelanggan:

5 orang memilih kopi

5 orang memilih teh

5 orang memilih cokelat

Padahal secara probabilitas, kopi dan teh lebih disukai (masing-masing 40%) dibandingkan cokelat (20%). Namun, hasil simulasi ini menunjukkan bahwa dalam sampel acak, jumlah pilihan tiap minuman bisa saja sama meskipun probabilitasnya berbeda. Ini adalah karakteristik dari distribusi multinomial, yang memodelkan hasil percobaan acak dengan lebih dari dua kategori.

Distribusi Poisson

Distribusi Poisson digunakan untuk memodelkan jumlah kejadian yang terjadi dalam suatu interval waktu atau ruang, ketika kejadian-kejadian tersebut: - Terjadi secara acak, - Terpisah satu sama lain (tidak berkelompok), - Memiliki rata-rata kejadian yang tetap (dilambangkan dengan \(\lambda\)).

Nilai \(\lambda\) menyatakan rata-rata jumlah kejadian per unit (misalnya per jam, per kilometer, atau per meter persegi). Distribusi ini cocok digunakan ketika kita ingin menghitung peluang banyaknya kejadian dalam rentang yang ditentukan.

Fungsi probabilitas dari distribusi Poisson:

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

Keterangan: - \(X\): variabel acak yang menyatakan jumlah kejadian - \(\lambda\): rata-rata jumlah kejadian dalam interval tertentu - \(e\): bilangan eksponensial (sekitar 2.718) - \(k\): jumlah kejadian yang diamati (\(k = 0, 1, 2, \dots\))

Contoh Kasus

Jumlah telepon masuk ke call center dalam satu jam. Rata-rata telepon yang masuk adalah 4 per jam. Simulasikan 10 observasi dari distribusi Poisson dengan \(\lambda = 4\).

set.seed(123)
poisson_sample <- rpois(10, lambda = 4)
poisson_sample
##  [1] 3 6 3 6 7 1 4 7 4 4

Kesimpulan

Hasil di atas menunjukkan jumlah telepon masuk ke call center selama 10 jam. Nilai-nilai tersebut dapat diartikan sebagai berikut:

Jam 1: 3 telepon

Jam 2: 6 telepon

Jam 3: 3 telepon

Jam 4: 6 telepon

Jam 5: 7 telepon

Jam 6: 1 telepon

Jam 7: 4 telepon

Jam 8: 7 telepon

Jam 9: 4 telepon

Jam 10: 4 telepon

Meskipun rata-rata jumlah telepon adalah 4 per jam, jumlah aktual pada tiap jam bervariasi — kadang di bawah, kadang di atas 4. Ini adalah ciri khas dari distribusi Poisson, yang cocok untuk memodelkan jumlah kejadian dalam waktu tertentu dengan rata-rata tetap, namun dengan kejadian yang bersifat acak.

Desain Sampling dalam Analisis Data Kategori

Desain sampling merupakan tahap penting dalam proses pengumpulan data, karena menentukan bagaimana data diambil dari populasi dan berdampak langsung pada validitas serta reliabilitas hasil penelitian. Dalam analisis data kategori, desain sampling juga memengaruhi cara hubungan antar variabel dianalisis dan disimpulkan secara tepat.

Secara umum, terdapat dua pendekatan utama dalam desain sampling untuk data kategori, yaitu sampling prospektif dan sampling retrospektif. Kedua pendekatan ini memiliki karakteristik serta metode pengambilan sampel yang berbeda, dan sering digunakan dalam berbagai jenis penelitian, baik eksperimen, studi kohort, maupun studi kasus-kontrol, tergantung pada tujuan dan jenis data yang dikumpulkan.

Desain Prospektif

Desain prospektif dilakukan dengan mengidentifikasi kelompok berdasarkan paparan terhadap suatu faktor (exposure), kemudian mengikuti mereka seiring waktu untuk mengamati apakah mereka mengalami kejadian tertentu (outcome). Pendekatan ini sering digunakan dalam penelitian eksperimental maupun observasional, seperti studi kohort, karena memungkinkan peneliti untuk menganalisis hubungan sebab-akibat secara lebih akurat.

Eksperimen

Dalam studi eksperimental, subjek penelitian dibagi secara acak ke dalam kelompok perlakuan dan kontrol. Tujuan utamanya adalah untuk mengamati efek dari suatu intervensi secara langsung. Beberapa teknik sampling yang umum digunakan dalam eksperimen meliputi:

  • Simple Random Sampling (SRS)
    Setiap individu dalam populasi memiliki peluang yang sama untuk dipilih sebagai sampel.

  • Stratified Random Sampling
    Populasi dibagi ke dalam strata berdasarkan karakteristik tertentu, lalu sampel diambil secara acak dari masing-masing strata untuk memastikan keterwakilan yang seimbang.

  • Cluster Sampling
    Populasi dibagi ke dalam beberapa klaster, lalu beberapa klaster dipilih secara acak dan seluruh individu di dalam klaster tersebut dianalisis.

Studi Kohort

Studi kohort merupakan jenis penelitian observasional di mana sekelompok individu dengan karakteristik tertentu (biasanya berdasarkan paparan terhadap suatu faktor) diikuti dari waktu ke waktu untuk mengamati perkembangan kejadian tertentu. Teknik sampling yang umum digunakan dalam studi kohort antara lain:

  • Census Sampling
    Seluruh anggota populasi yang memenuhi kriteria penelitian dilibatkan.

  • Systematic Sampling
    Subjek dipilih dari populasi berdasarkan interval tertentu (misalnya setiap orang ke-10 dari daftar).

  • Matched Sampling
    Setiap individu dalam kelompok kohort dicocokkan dengan individu serupa dari kelompok pembanding berdasarkan variabel-variabel tertentu untuk mengurangi bias.

Contoh Kasus

Peneliti ingin mengetahui apakah konsumsi kopi berlebihan dapat meningkatkan risiko insomnia. Untuk itu, mereka mengidentifikasi dua kelompok:

  • Kelompok pertama terdiri dari individu yang mengonsumsi lebih dari 3 cangkir kopi per hari.
  • Kelompok kedua terdiri dari individu yang tidak mengonsumsi kopi dalam jumlah tersebut.

Kedua kelompok tersebut diikuti selama 30 hari untuk melihat berapa banyak individu dalam masing-masing kelompok yang mengalami insomnia. Karena data dikumpulkan ke depan dari waktu pengamatan awal, maka studi ini termasuk dalam desain prospektif.

Desain Retrospektif

Desain retrospektif digunakan ketika data dikumpulkan dari kejadian yang telah terjadi di masa lalu. Pendekatan ini sering digunakan dalam studi observasional yang bertujuan untuk menelusuri hubungan antara faktor risiko dan outcome berdasarkan informasi historis. Karena data dikumpulkan setelah outcome terjadi, penelitian ini cenderung lebih cepat dan hemat biaya dibandingkan desain prospektif.

Studi Kasus-Kontrol

Dalam studi kasus-kontrol, peneliti memilih dua kelompok: satu kelompok dengan kondisi tertentu (kasus) dan satu kelompok tanpa kondisi tersebut (kontrol). Peneliti kemudian menelusuri riwayat paparan atau karakteristik masa lalu untuk mencari hubungan antara paparan dan outcome. Teknik sampling yang sering digunakan dalam studi ini meliputi:

  • Purposive Sampling
    Pemilihan sampel berdasarkan karakteristik yang dianggap relevan dengan tujuan penelitian.

  • Snowball Sampling
    Sampel awal membantu merekrut individu lain yang memiliki karakteristik serupa, biasanya digunakan saat populasi target sulit dijangkau.

  • Incidence Density Sampling
    Sampel kasus dan kontrol dipilih dari populasi yang sama, disesuaikan berdasarkan waktu kemunculan kasus untuk mengontrol potensi bias waktu.

Studi Kohort Retrospektif

Studi kohort retrospektif dilakukan dengan menggunakan data historis untuk mengelompokkan individu berdasarkan paparan masa lalu, lalu menganalisis apakah mereka mengalami outcome tertentu. Berbeda dari studi kohort prospektif, pengamatan dilakukan ke belakang. Beberapa metode sampling yang umum digunakan:

  • Convenience Sampling
    Sampel dipilih berdasarkan ketersediaan data, bukan representasi populasi secara acak.

  • Quota Sampling
    Sampel dipilih untuk mencerminkan proporsi tertentu dari populasi, berdasarkan karakteristik tertentu.

  • Case-Based Sampling
    Sampel ditentukan berdasarkan kejadian yang telah terjadi, misalnya berdasarkan data kasus di rumah sakit atau catatan medis.

Contoh Kasus

Seorang peneliti ingin menyelidiki hubungan antara paparan pestisida dan kejadian kanker. Ia menggunakan data rekam medis dan riwayat pekerjaan petani yang telah terdiagnosis kanker (kasus) dan membandingkannya dengan kelompok petani sehat (kontrol). Karena baik paparan maupun outcome telah terjadi sebelum penelitian dimulai, studi ini merupakan retrospektif.

Tabel Perbandingan Desain Sampling

Tabel berikut merangkum perbandingan berbagai jenis desain sampling yang umum digunakan dalam analisis data kategori. Tabel ini mencakup pendekatan yang digunakan, metode sampling, serta kelebihan dan kekurangannya:

Jenis Studi Pendekatan Metode Sampling Keuntungan Kelemahan
Eksperimen Prospektif SRS, Stratified, Cluster Kontrol tinggi terhadap variabel, hasil dapat dianalisis dengan akurat Biaya tinggi, etika dan validitas perlu diperhatikan
Studi Kohort Prospektif Census, Systematic, Matched Dapat mengamati perkembangan kejadian dalam jangka panjang Membutuhkan waktu lama, risiko kehilangan partisipan
Studi Kasus-Kontrol Retrospektif Purposive, Snowball, Incidence Density Mudah dan cepat dilakukan, efisien untuk penyakit langka Sulit mengontrol variabel pengganggu, rentan bias recall
Studi Kohort Retrospektif Retrospektif Convenience, Quota, Case-Based Memanfaatkan data historis, lebih murah daripada studi prospektif Kualitas data bervariasi, bias data historis, risiko missing data

Kesimpulan

Desain sampling dalam analisis data kategori sangat dipengaruhi oleh pendekatan yang digunakan, apakah itu prospektif atau retrospektif. Pemilihan metode sampling yang sesuai tidak hanya bergantung pada jenis studi, tetapi juga tujuan penelitian, ketersediaan data, dan sumber daya yang ada. Studi prospektif seperti eksperimen dan kohort lebih kuat secara validitas kausalitas, namun membutuhkan waktu dan biaya lebih besar. Di sisi lain, studi retrospektif seperti kasus-kontrol dan kohort retrospektif lebih efisien dari segi waktu dan biaya, tetapi rentan terhadap bias data historis. Oleh karena itu, penting untuk memilih desain dan metode sampling yang tepat agar menghasilkan analisis yang valid dan dapat dipertanggungjawabkan.

Tabel Kontingensi 2x2

Tabel kontingensi 2×2 adalah jenis tabel yang paling dasar untuk menganalisis hubungan antara dua variabel kategori. Tabel ini sering digunakan dalam analisis statistik untuk menguji apakah ada hubungan atau asosiasi antara dua variabel, misalnya hubungan antara jenis pengobatan dan tingkat kesembuhan, atau antara kebiasaan merokok dengan risiko terkena kanker paru-paru. Struktur tabel kontingensi 2×2 sangat sederhana dan memuat data yang membandingkan dua variabel kategori, masing-masing dengan dua kategori.

Struktur umum tabel ini mencerminkan jumlah kejadian dan non-kejadian pada dua kelompok (misalnya: kelompok yang terpapar dan tidak terpapar).

Struktur Umum

Struktur tabel kontingensi 2x2 dapat ditulis sebagai berikut:

# Buat data frame
df <- data.frame(
  Kejadian_Plus = c("n11", "n21", "n.1"),
  Kejadian_Minus = c("n12", "n22", "n.2"),
  Total = c("n1.", "n2.", "n"),
  row.names = c("Paparan (+)", "Paparan (-)", "Total")
)

# Tampilkan tabel
knitr::kable(
  df,
  caption = "Struktur Umum Tabel Kontingensi 2x2"
)
Struktur Umum Tabel Kontingensi 2x2
Kejadian_Plus Kejadian_Minus Total
Paparan (+) n11 n12 n1.
Paparan (-) n21 n22 n2.
Total n.1 n.2 n
  • \(n_{11}\): Jumlah kasus dalam kategori positif dari Grup 1
  • \(n_{12}\): Jumlah kasus dalam kategori negatif dari Grup 1
  • \(n_{21}\): Jumlah kasus dalam kategori positif dari Grup 2
  • \(n_{22}\): Jumlah kasus dalam kategori negatif dari Grup 2
  • \(n_{1.} = n_{11} + n_{12}\): Total observasi dalam Grup 1
  • \(n_{2.} = n_{21} + n_{22}\): Total observasi dalam Grup 2
  • \(n_{.1} = n_{11} + n_{21}\): Total observasi dalam Kategori Kejadian_Plus
  • \(n_{.2} = n_{12} + n_{22}\): Total observasi dalam Kategori Kejadian_Minus
  • \(n = n_{11} + n_{12} + n_{21} + n_{22}\): Total keseluruhan sampel

Contoh Kasus

Misalkan dilakukan studi terhadap 200 individu untuk mengetahui apakah vaksinasi menurunkan risiko infeksi penyakit X.

  • Dari 90 orang yang divaksinasi, 10 orang terinfeksi.
  • Dari 110 orang tanpa vaksinasi, 30 orang terinfeksi.

Tabel kontingensi-nya adalah sebagai berikut:

data <- matrix(c(10, 80, 30, 80), nrow = 2, byrow = TRUE)
colnames(data) <- c("Terinfeksi", "Tidak Terinfeksi")
rownames(data) <- c("Vaksinasi", "Tanpa Vaksin")
knitr::kable(data, caption = "Tabel Kontingensi: Status Vaksinasi dan Infeksi")
Tabel Kontingensi: Status Vaksinasi dan Infeksi
Terinfeksi Tidak Terinfeksi
Vaksinasi 10 80
Tanpa Vaksin 30 80

Distribusi Peluang dalam Tabel Kontingensi 2 × 2

Tabel kontingensi 2×2 sering digunakan untuk mengidentifikasi hubungan antara dua variabel kategori. Distribusi peluang dalam tabel ini bisa dianalisis dengan tiga jenis peluang:

  • Peluang Bersama
  • Peluang Marginal
  • Peluang Bersyarat

Sebagai contoh, kita akan menganalisis hubungan antara kebiasaan konsumsi soda dan kejadian obesitas berdasarkan data survei terhadap 800 responden.

Peluang Bersama

Peluang bersama mengacu pada peluang terjadinya dua peristiwa secara bersamaan. Dalam tabel kontingensi 2×2, peluang bersama dihitung dengan membagi jumlah kasus di suatu sel dengan total seluruh observasi.

Secara matematis, peluang bersama untuk dua peristiwa A dan B dapat dihitung dengan rumus:

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

Contoh Kasus:

Kebiasaan Obesitas (+) Obesitas (-) Total
Minum Soda 280 120 400
Tidak Minum Soda 70 330 400
Total 350 450 800 
# Membuat tabel kontingensi
tabel <- matrix(c(280, 120, 70, 330), 
                nrow = 2, byrow = TRUE,
                dimnames = list(
                  "Kebiasaan" = c("Minum Soda", "Tidak Minum Soda"),
                  "Obesitas" = c("Ya", "Tidak")
                ))
total <- sum(tabel)
tabel
##                   Obesitas
## Kebiasaan           Ya Tidak
##   Minum Soda       280   120
##   Tidak Minum Soda  70   330

Hitung Peluang Bersama untuk tiap sel:

peluang_bersama <- tabel / total
peluang_bersama
##                   Obesitas
## Kebiasaan              Ya  Tidak
##   Minum Soda       0.3500 0.1500
##   Tidak Minum Soda 0.0875 0.4125

Interpretasi Peluang Bersama:

  • \(P(\text{Minum Soda} \cap \text{Obesitas}) = \frac{280}{800} = 0.35\)
  • \(P(\text{Minum Soda} \cap \text{Tidak Obesitas}) = \frac{120}{800} = 0.15\)
  • \(P(\text{Tidak Minum} \cap \text{Obesitas}) = \frac{70}{800} = 0.0875\)
  • \(P(\text{Tidak Minum} \cap \text{Tidak Obesitas}) = \frac{330}{800} = 0.4125\)

Peluang Marginal

Peluang marginal adalah probabilitas terjadinya suatu peristiwa tunggal, tanpa memperhitungkan peristiwa lainnya.

Dalam konteks tabel kontingensi 2×2, peluang marginal diperoleh dengan:

  • Menjumlahkan frekuensi baris tertentu (untuk variabel baris), lalu membaginya dengan total keseluruhan.
  • Atau menjumlahkan frekuensi kolom tertentu (untuk variabel kolom), lalu membaginya dengan total keseluruhan.

Secara matematis:

  • Untuk baris: \(P(A_i) = \frac{n_{i.}}{n}\)
  • Untuk kolom: \(P(B_j) = \frac{n_{.j}}{n}\)

Keterangan: - \(n_{i.}\): jumlah frekuensi seluruh kolom pada baris ke-\(i\) - \(n_{.j}\): jumlah frekuensi seluruh baris pada kolom ke-\(j\) - \(n\): total seluruh pengamatan

# Peluang marginal berdasarkan baris (kebiasaan)
peluang_marginal_baris <- rowSums(tabel) / total
peluang_marginal_baris
##       Minum Soda Tidak Minum Soda 
##              0.5              0.5
# Peluang marginal berdasarkan kolom (obesitas)
peluang_marginal_kolom <- colSums(tabel) / total
peluang_marginal_kolom
##     Ya  Tidak 
## 0.4375 0.5625

Interpretasi:

  • \(P(\text{Minum Soda}) = \frac{400}{800} = 0.5\)
  • \(P(\text{Tidak Minum Soda}) = \frac{400}{800} = 0.5\)
  • \(P(\text{Obesitas}) = \frac{350}{800} = 0.4375\)
  • \(P(\text{Tidak Obesitas}) = \frac{450}{800} = 0.5625\)

Peluang Bersyarat

Peluang bersyarat mengukur kemungkinan suatu peristiwa terjadi, dengan syarat peristiwa lainnya sudah terjadi. Dalam tabel kontingensi 2×2, peluang bersyarat dihitung dengan membagi peluang bersama dengan peluang marginal dari peristiwa yang bersyarat.

Secara matematis :

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

# Peluang bersyarat: Obesitas jika Minum Soda
p_obesitas_given_soda <- tabel["Minum Soda", "Ya"] / sum(tabel["Minum Soda", ])
p_obesitas_given_soda
## [1] 0.7
# Peluang bersyarat: Tidak Obesitas jika Minum Soda
p_tidak_obesitas_given_soda <- tabel["Minum Soda", "Tidak"] / sum(tabel["Minum Soda", ])
p_tidak_obesitas_given_soda
## [1] 0.3
# Peluang bersyarat: Obesitas jika Tidak Minum Soda
p_obesitas_given_no_soda <- tabel["Tidak Minum Soda", "Ya"] / sum(tabel["Tidak Minum Soda", ])
p_obesitas_given_no_soda
## [1] 0.175
# Peluang bersyarat: Tidak Obesitas jika Tidak Minum Soda
p_tidak_obesitas_given_no_soda <- tabel["Tidak Minum Soda", "Tidak"] / sum(tabel["Tidak Minum Soda", ])
p_tidak_obesitas_given_no_soda
## [1] 0.825

Interpretasi:

  • \(P(\text{Obesitas} \mid \text{Minum Soda}) = \frac{280}{400} = 0.70\)
  • \(P(\text{Tidak Obesitas} \mid \text{Minum Soda}) = \frac{120}{400} = 0.30\)
  • \(P(\text{Obesitas} \mid \text{Tidak Minum Soda}) = \frac{70}{400} = 0.175\)
  • \(P(\text{Tidak Obesitas} \mid \text{Tidak Minum Soda}) = \frac{330}{400} = 0.825\)

Kesimpulan

Dari analisis peluang di atas:

  • Konsumsi soda tampaknya berhubungan dengan tingkat obesitas yang lebih tinggi.
  • Peluang bersyarat membantu kita memahami efek suatu faktor terhadap variabel lain.
Tabel Kontingensi: Konsumsi Soda dan Obesitas
Ya Tidak
Minum Soda 280 120
Tidak Minum Soda 70 330

Ukuran Asosiasi dalam Data Kategori 2 × 2

Analisis asosiasi dalam tabel kontingensi 2 × 2 merupakan pendekatan statistik yang digunakan untuk mengevaluasi hubungan antara dua variabel kategori. Tabel ini menampilkan frekuensi observasi berdasarkan kombinasi dua variabel, masing-masing dengan dua kategori. Karena kesederhanaannya, tabel ini banyak digunakan untuk memahami keterkaitan atau pengaruh antar variabel.

Tabel 2 × 2 sangat umum digunakan dalam:

  • Epidemiologi: Hubungan merokok dengan kanker paru.
  • Klinis: Efektivitas vaksin.
  • Sosial: Hubungan pendidikan dengan pengangguran.

Contoh struktur tabel:

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

Risk Difference (RD)

Risk Difference (RD) atau selisih risiko adalah ukuran absolut yang mengukur perbedaan proporsi kejadian antar dua kelompok. RD menghitung seberapa besar risiko lebih (atau kurang) dalam kelompok yang terpapar dibandingkan kelompok yang tidak terpapar.

\[ RD = \frac{n_{11}}{n_{1.}} - \frac{n_{21}}{n_{2.}} \]

  • Jika \(RD > 0\), maka risiko kejadian lebih tinggi di kelompok pertama.
  • Jika \(RD < 0\), maka risiko kejadian lebih rendah di kelompok pertama.
  • Jika \(RD = 0\), maka tidak ada perbedaan risiko.

Contoh Kasus:

Sebuah penelitian ingin mengetahui apakah konsumsi makanan cepat saji berkaitan dengan kejadian kolesterol tinggi. Dari 200 responden:

  • 90 orang konsumsi fast food secara rutin, dan 40 di antaranya punya kolesterol tinggi.
  • 110 orang tidak konsumsi fast food secara rutin, dan 22 di antaranya punya kolesterol tinggi.

Tabel Kontingensi:

Kolesterol Tinggi Tidak Tinggi Total
Fast Food Rutin 40 50 90
Tidak Fast Food 22 88 110
Total 62 138 200 
RD <- function(n11, n12, n21, n22) {
  (n11 / (n11 + n12)) - (n21 / (n21 + n22))
}
RD(40, 50, 22, 88)
## [1] 0.2444444

Interpretasi:

Risk Difference menunjukkan bahwa terdapat perbedaan risiko antara dua kelompok. Nilai positif berarti kelompok fast food rutin lebih berisiko mengalami kolesterol tinggi.

Relative Risk (RR)

Relative Risk (RR) atau risiko relatif menunjukkan rasio antara probabilitas kejadian pada kelompok terpapar dibandingkan dengan yang tidak terpapar. RR menunjukkan “berapa kali lebih mungkin” kejadian terjadi pada kelompok terpapar.

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

  • \(RR > 1\): kejadian lebih sering terjadi di kelompok pertama.
  • \(RR < 1\): kejadian lebih jarang terjadi di kelompok pertama.
  • \(RR = 1\): tidak ada perbedaan relatif antara kedua kelompok.
RR <- function(n11, n12, n21, n22) {
  (n11 / (n11 + n12)) / (n21 / (n21 + n22))
}
RR(40, 50, 22, 88)
## [1] 2.222222

Interpretasi:

Relative Risk mengukur seberapa besar kemungkinan mengalami kolesterol tinggi jika rutin makan fast food, dibandingkan dengan yang tidak.

Odds Ratio (OR)

Odds Ratio (OR) adalah perbandingan antara peluang (odds) kejadian di satu kelompok dengan peluang di kelompok lainnya. OR sering digunakan dalam studi kasus-kontrol karena tidak membutuhkan total populasi.

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

  • \(OR > 1\): odds lebih besar di kelompok pertama.
  • \(OR < 1\): odds lebih kecil di kelompok pertama.
  • \(OR = 1\): tidak ada perbedaan odds antara kelompok.
OR <- function(n11, n12, n21, n22) {
  (n11 * n22) / (n12 * n21)
}
OR(40, 50, 22, 88)
## [1] 3.2

Interpretasi:

Odds Ratio memperlihatkan perbandingan peluang antara kelompok fast food dan tidak fast food untuk mengalami kolesterol tinggi.

Tabel Perbandingan RD, RR, dan OR

Ukuran Definisi Studi yang Cocok Interpretasi
Risk Difference Selisih proporsi kejadian antar dua kelompok Studi kohort, eksperimen Mengukur dampak absolut suatu faktor
Relative Risk Rasio probabilitas kejadian antar dua kelompok Studi kohort, klinis Mengukur risiko relatif (berapa kali)
Odds Ratio Rasio peluang kejadian antar dua kelompok Studi kasus-kontrol, observasional Mengukur odds antar dua kelompok

Kesimpulan:
- RD cocok untuk mengevaluasi efek absolut dari suatu intervensi.
- RR umum digunakan untuk studi prospektif karena mempertimbangkan probabilitas.
- OR sering dipakai di studi observasional atau kasus-kontrol karena tidak butuh total populasi.

Inferensi Tabel Kontingensi Dua Arah

Tabel kontingensi dua arah menyajikan frekuensi gabungan dari dua variabel kategori dalam bentuk matriks. Tujuannya adalah untuk melihat apakah distribusi satu variabel dipengaruhi oleh variabel lainnya.

Inferensi dalam tabel ini terbagi menjadi dua bagian utama:

Dengan proses ini, kita dapat menilai apakah ada hubungan yang berarti antara dua variabel kategori berdasarkan data yang ada.

Menguji apakah terdapat perbedaan proporsi signifikan antara dua kelompok.

Contoh Kasus:

Penelitian terhadap dua metode pembelajaran: 65 dari 100 siswa lulus pada metode A, dan 55 dari 100 siswa lulus pada metode B.

Estimasi

Estimasi dalam statistik digunakan untuk memberikan informasi tentang parameter populasi berdasarkan data sampel. Estimasi ini dibagi menjadi dua jenis utama, yaitu estimasi titik dan estimasi interval.

Estimasi Titik

Estimasi titik bertujuan memberikan satu nilai terbaik sebagai perkiraan terhadap parameter populasi. Dalam konteks proporsi, estimasi titik dapat dihitung dengan:

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

Di mana: - \(\hat{p}\) adalah estimasi titik proporsi, - \(x\) adalah jumlah individu dalam kategori tertentu, - \(n\) adalah jumlah total individu dalam sampel.

Contoh Kasus:

Misalnya, dari 100 orang yang mengikuti program diet, 62 orang berhasil menurunkan berat badan. Maka estimasi titik proporsi keberhasilan adalah:

x <- 62
n <- 100
p_hat <- x/n
p_hat
## [1] 0.62

Artinya, diperkirakan 62% peserta berhasil menurunkan berat badan.

Estimasi Interval

Estimasi interval bertujuan memberikan rentang nilai yang dipercaya memuat nilai sebenarnya dari parameter populasi dengan tingkat kepercayaan tertentu. Rumus umum untuk proporsi:

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

Contoh Kasus: Gunakan data sebelumnya dengan tingkat kepercayaan 95% (\(Z_{\alpha/2} = 1.96\)):

z <- 1.96
se <- sqrt(p_hat * (1 - p_hat) / n)
lower <- p_hat - z * se
upper <- p_hat + z * se
c(lower, upper)
## [1] 0.5248643 0.7151357

Interpretasi

Kita 95% yakin bahwa proporsi keberhasilan program diet berada di antara batas bawah dan atas yang dihitung.

Uji Hipotesis

Uji hipotesis digunakan untuk menguji dugaan terhadap parameter populasi berdasarkan data sampel. Dalam konteks tabel kontingensi 2x2, fokus kita adalah pada uji proporsi.

Uji Proporsi

Uji proporsi dua kelompok digunakan untuk membandingkan apakah proporsi suatu kejadian pada dua kelompok berbeda secara signifikan.

Misalnya kita memiliki data sebagai berikut:

Kejadian (+) Tidak Kejadian (-) Total
Terpapar 50 30 80
Tidak Terpapar 30 50 80 
set.seed(123)
data <- matrix(c(50, 30, 30, 50), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Terpapar" = c("Ya", "Tidak"), "Kejadian" = c("Ya", "Tidak"))
data
##         Kejadian
## Terpapar Ya Tidak
##    Ya    50    30
##    Tidak 30    50

Gunakan fungsi prop.test() untuk menguji apakah ada perbedaan proporsi yang signifikan:

prop_test <- prop.test(x = c(data[1,1], data[2,1]),
                       n = c(sum(data[1,]), sum(data[2,])))
prop_test
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(data[1, 1], data[2, 1]) out of c(sum(data[1, ]), sum(data[2, ]))
## X-squared = 9.025, df = 1, p-value = 0.002663
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.08747151 0.41252849
## sample estimates:
## prop 1 prop 2 
##  0.625  0.375

Interpretasi:

  • Jika nilai p (p-value) < 0.05, maka kita menolak hipotesis nol dan menyimpulkan bahwa terdapat perbedaan proporsi yang signifikan antara kelompok terpapar dan tidak terpapar.
  • Jika p-value > 0.05, maka tidak ada bukti yang cukup untuk menyatakan bahwa proporsinya berbeda secara signifikan.

Kesimpulan

Karena p-value (0.007699) < 0.05, maka kita menolak hipotesis nol dan menyimpulkan bahwa terdapat perbedaan proporsi yang signifikan antara kelompok terpapar dan tidak terpapar.

Uji Asosiasi

Uji asosiasi dalam tabel kontingensi 2×2 bertujuan untuk menilai apakah terdapat hubungan atau ketergantungan antara dua variabel kategori.

Hipotesis Uji dalam Tabel Kontingensi 2×2

Untuk setiap uji asosiasi dalam tabel kontingensi 2×2, hipotesis yang diuji adalah sebagai berikut:

  1. Hipotesis Nol (H₀) Hipotesis nol menyatakan bahwa tidak ada asosiasi antara dua variabel. Dalam hal ini, distribusi frekuensi yang diamati pada tabel kontingensi 2×2 tidak berbeda secara signifikan dari distribusi yang diharapkan jika kedua variabel tersebut independen satu sama lain.

\[ H_0: \text{Tidak ada asosiasi antara kedua variabel.} \]

  1. Hipotesis Alternatif (H₁) Hipotesis alternatif menyatakan bahwa terdapat asosiasi antara dua variabel. Artinya, distribusi frekuensi yang diamati pada tabel kontingensi 2×2 berbeda secara signifikan dari distribusi yang diharapkan jika kedua variabel tersebut tergantung satu sama lain.

\[ H_1: \text{Terdapat asosiasi antara kedua variabel.} \]

  1. Interpretasi Hasil
  • Jika p-value < 0.05, kita menolak hipotesis nol dan menyimpulkan ada asosiasi antara kedua variabel.
  • Jika p-value ≥ 0.05, kita gagal menolak hipotesis nol dan menyimpulkan bahwa tidak ada asosiasi yang signifikan antara kedua variabel.

Risk Difference (RD)

Mengukur selisih antara risiko kejadian suatu peristiwa pada dua kelompok. RD dihitung dengan cara mengurangi risiko pada kelompok kontrol dari risiko pada kelompok eksperimen. Ini memberikan gambaran seberapa besar perbedaan risiko absolut antara dua kelompok, misalnya dalam penelitian klinis untuk melihat efek pengobatan.

\[ 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)} \]

Contoh Kasus:

Sebuah studi membandingkan efek vaksin terhadap kejadian infeksi. Dari 150 orang yang divaksin, 15 mengalami infeksi. Dari 130 orang yang tidak divaksin, 39 mengalami infeksi.

n11 <- 15
n12 <- 135
n21 <- 39
n22 <- 91

n1. <- n11 + n12
n2. <- n21 + n22

p1 <- n11 / n1.
p2 <- n21 / n2.
rd <- p1 - p2
se_rd <- sqrt((p1 * (1 - p1) / n1.) + (p2 * (1 - p2) / n2.))
z_rd <- rd / se_rd

list(RD = rd, SE_RD = se_rd, Z_RD = z_rd)
## $RD
## [1] -0.2
## 
## $SE_RD
## [1] 0.04706787
## 
## $Z_RD
## [1] -4.249183

Kesimpulan

Risk Difference (RD) menunjukkan selisih risiko infeksi antara kelompok yang divaksin dan tidak divaksin. Jika hasil RD bernilai negatif, berarti vaksin mengurangi risiko infeksi dibandingkan kelompok yang tidak divaksin. Misalnya jika diperoleh RD = -0.17, maka terdapat pengurangan risiko infeksi sebesar 17% pada kelompok yang divaksin dibandingkan yang tidak. Jika nilai Z signifikan (misalnya |Z| > 1.96 pada α = 0.05), maka perbedaan ini secara statistik signifikan.

Relative Risk (RR)

Mengukur perbandingan risiko antara dua kelompok. RR dihitung dengan membandingkan risiko kejadian pada kelompok eksperimen dengan risiko kejadian pada kelompok kontrol. Jika RR > 1, menunjukkan bahwa kelompok eksperimen memiliki risiko yang lebih tinggi, sedangkan jika RR < 1, menunjukkan risiko yang lebih rendah dibandingkan kontrol.

\[ RR = \frac{\frac{n_{11}}{n_{1.}}}{\frac{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)} \]

Contoh Kasus:

Dalam uji coba obat baru, dari 100 pasien yang diberi obat, 20 sembuh. Dari 120 pasien yang diberi plasebo, hanya 15 sembuh.

n11 <- 20
n12 <- 80
n21 <- 15
n22 <- 105

n1. <- n11 + n12
n2. <- n21 + n22

rr <- (n11 / n1.) / (n21 / n2.)
se_ln_rr <- sqrt((1 / n11) - (1 / n1.) + (1 / n21) - (1 / n2.))
z_rr <- log(rr) / se_ln_rr

list(RR = rr, SE_Ln_RR = se_ln_rr, Z_RR = z_rr)
## $RR
## [1] 1.6
## 
## $SE_Ln_RR
## [1] 0.3135815
## 
## $Z_RR
## [1] 1.498825

Kesimpulan

Relative Risk (RR) membandingkan kemungkinan sembuh antara pasien yang menerima obat dengan pasien yang menerima plasebo. Jika RR > 1, maka kemungkinan sembuh lebih besar pada kelompok yang diberi obat. Misalnya jika RR = 1.6, berarti pasien yang diberi obat 1.6 kali lebih mungkin sembuh dibandingkan yang diberi plasebo. Jika Z signifikan, maka perbedaan risiko ini tidak terjadi secara kebetulan.

Odds Ratio (OR)

Mengukur perbandingan odds (peluang) kejadian suatu peristiwa antara dua kelompok. Odds adalah perbandingan antara kemungkinan terjadinya dan tidak terjadinya peristiwa. OR sering digunakan dalam studi kasus kontrol untuk melihat kekuatan asosiasi antara faktor risiko dan kejadian penyakit.

\[ 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)} \]

Contoh Kasus:

Sebuah penelitian ingin mengetahui hubungan antara kebiasaan merokok dan kejadian kanker paru. Dari 60 perokok, 18 mengidap kanker paru, sedangkan dari 80 non-perokok, 10 mengalami kanker.

n11 <- 18
n12 <- 42
n21 <- 10
n22 <- 70

n1. <- n11 + n12
n2. <- n21 + n22

or <- (n11 * n22) / (n12 * n21)
se_ln_or <- sqrt((1 / n11) + (1 / n12) + (1 / n21) + (1 / n22))
z_or <- log(or) / se_ln_or

list(OR = or, SE_Ln_OR = se_ln_or, Z_OR = z_or)
## $OR
## [1] 3
## 
## $SE_Ln_OR
## [1] 0.4400577
## 
## $Z_OR
## [1] 2.496519

Kesimpulan

Odds Ratio (OR) menunjukkan perbandingan peluang terjadinya kanker paru antara perokok dan non-perokok. Jika OR > 1, maka perokok memiliki peluang lebih tinggi untuk terkena kanker paru dibandingkan non-perokok. Misalnya jika OR = 2.8, maka perokok memiliki odds terkena kanker paru hampir 3 kali lebih tinggi dibandingkan non-perokok. Jika nilai Z signifikan, maka hubungan antara merokok dan kanker paru secara statistik penting.

Uji Independensi

Uji independensi digunakan untuk mengetahui apakah ada hubungan antara dua variabel kategorik. Salah satu metode paling umum adalah uji Chi-Square. Selain itu, terdapat juga uji alternatif seperti Likelihood Ratio (G²) dan uji Fisher yang digunakan ketika asumsi Chi-Square tidak terpenuhi.

Uji Chi-Square

Uji Chi-Square (χ²) digunakan untuk mengetahui apakah dua variabel kategorik bersifat independen atau tidak. Rumus dari statistik Chi-Square adalah:

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

Di mana: - \(O\): Frekuensi observasi (data yang sebenarnya terjadi) - \(E\): Frekuensi harapan (data yang diharapkan jika variabel bebas), dihitung dengan:

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

Contoh Kasus:

Sebuah survei dilakukan untuk melihat hubungan antara penggunaan masker dan kejadian flu.

# Data Kontingensi
data <- matrix(c(30, 10, 15, 45), nrow = 2, byrow = TRUE)
dimnames(data) <- list("Masker" = c("Ya", "Tidak"),
                       "Flu" = c("Ya", "Tidak"))
data
##        Flu
## Masker  Ya Tidak
##   Ya    30    10
##   Tidak 15    45
# Uji Chi-Square
chisq_test <- chisq.test(data)
chisq_test
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data
## X-squared = 22.264, df = 1, p-value = 2.376e-06

Interpretasi

Nilai p-value yang dihasilkan lebih kecil dari 0.05 berarti kita menolak H0, artinya terdapat hubungan antara penggunaan masker dan kejadian flu.

Kesimpulan

Karena p-value (2.376e-06) > 0.05, maka menerima H0, artinya tidak terdapat hubungan natara penggunaan masker dan kejadian flu dalam kasus ini.

Partisi Chi-Square

Partisi Chi-Square digunakan untuk mengidentifikasi lebih lanjut pada tabel yang kompleks (lebih dari 2x2). Dengan membagi tabel menjadi beberapa bagian (partisi), kita dapat melihat kontribusi masing-masing kategori terhadap hubungan yang terjadi.

Contoh Kasus:

Data mengenai jenis transportasi yang digunakan berdasarkan jenis kelamin.

# Tabel Kontingensi
data_matrix <- matrix(c(120, 90, 150, 100, 70, 180), nrow = 2, byrow = TRUE)
colnames(data_matrix) <- c("Mobil", "Motor", "Transportasi Umum")
rownames(data_matrix) <- c("Perempuan", "Laki-laki")
data_matrix
##           Mobil Motor Transportasi Umum
## Perempuan   120    90               150
## Laki-laki   100    70               180
# Uji Chi-Square keseluruhan
chi_test <- chisq.test(data_matrix)
chi_test
## 
##  Pearson's Chi-squared test
## 
## data:  data_matrix
## X-squared = 6.906, df = 2, p-value = 0.03165
# Partisi 1: Mobil vs Motor
data_part1 <- matrix(c(120, 90, 100, 70), nrow = 2, byrow = TRUE)
chisq.test(data_part1)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data_part1
## X-squared = 0.050832, df = 1, p-value = 0.8216
# Partisi 2: Gabungan Mobil+Motor vs Transportasi Umum
data_part2 <- matrix(c(210, 150, 180, 180), nrow = 2, byrow = TRUE)
chisq.test(data_part2)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data_part2
## X-squared = 4.7049, df = 1, p-value = 0.03008

Interpretasi:

Dengan mempartisi tabel, kita bisa mengetahui bahwa perbedaan signifikan berasal dari pilihan antara Mobil dan Motor, atau antara gabungan kendaraan pribadi dan transportasi umum.

  • Uji Chi-Square Tanpa Koreksi (Pearson’s Chi-squared test)

Interpretasi:

Nilai statistik uji: X-squared = 6.906

Derajat kebebasan: df = 2

Nilai-p: p-value = 0.03165

Karena nilai-p < 0.05, maka kita menolak H₀. Ini berarti terdapat hubungan yang signifikan secara statistik antara jenis kelamin dan pemilihan jenis transportasi.

  • Uji Chi-Square dengan Koreksi Yates (2x2)

Interpretasi: Nilai statistik uji: X-squared = 0.05083

df = 1

p-value = 0.8216

Karena nilai-p > 0.05, maka kita gagal menolak H₀. Artinya, tidak ada hubungan signifikan antara jenis kelamin dan pemilihan antara Mobil dan Motor.

  • Kasus 2: Mobil vs Transportasi Umum

Interpretasi: Nilai statistik uji: X-squared = 4.7049

df = 1

p-value = 0.03008

Karena nilai-p < 0.05, maka kita menolak H₀. Artinya, terdapat hubungan signifikan antara jenis kelamin dan pemilihan antara Mobil dan Transportasi Umum.

Kesimpulan

Secara keseluruhan, terdapat hubungan yang signifikan antara jenis kelamin dan pilihan transportasi.

Namun, hubungan ini tidak terlihat pada pilihan Mobil vs Motor, tetapi terlihat signifikan antara Mobil vs Transportasi Umum.

Uji Likelihood Ratio (G²)

Uji G² merupakan alternatif dari Chi-Square. G² menggunakan logaritma untuk mengukur perbedaan antara nilai observasi dan ekspektasi.

\[ G^2 = 2 \sum O \cdot \ln\left(\frac{O}{E}\right) \]

Contoh Kasus:

Penelitian dilakukan untuk melihat apakah merokok berkaitan dengan hipertensi.

# Data
data_matrix <- matrix(c(60, 40, 20, 80), nrow = 2, byrow = TRUE)
colnames(data_matrix) <- c("Hipertensi (+)", "Normal (-)")
rownames(data_matrix) <- c("Perokok", "Bukan Perokok")

# Ekspektasi dan G²
data_expected <- chisq.test(data_matrix)$expected
G2 <- 2 * sum(data_matrix * log(data_matrix / data_expected))
G2
## [1] 34.52185
# Nilai kritis untuk alpha = 0.05, df = 1
critical_value <- qchisq(0.95, df = 1)
critical_value
## [1] 3.841459
# Keputusan
if (G2 > critical_value) {
  "Tolak H0: Terdapat hubungan signifikan antara merokok dan hipertensi."
} else {
  "Gagal tolak H0: Tidak terdapat hubungan signifikan."
}
## [1] "Tolak H0: Terdapat hubungan signifikan antara merokok dan hipertensi."

Interpretasi:

  • Jika nilai G² lebih besar dari nilai kritis, maka kita menolak hipotesis nol dan menyimpulkan bahwa variabel tidak independen.

Kesimpulan

Karena G² = 34.52 > 3.84, maka keputusan uji adalah: Tolak H. Artinya, terdapat hubungan yang signifikan secara statistik antara status merokok dan kejadian hipertensi.

Orang yang merokok memiliki kecenderungan yang berbeda secara signifikan terhadap kemungkinan mengalami hipertensi dibandingkan dengan yang tidak merokok. Dengan kata lain, merokok berkorelasi dengan risiko hipertensi dalam data ini.

Uji Fisher’s Exact dan Distribusi Hipergeometrik

Jika ukuran sampel kecil dan tabel berukuran 2x2, uji Fisher’s Exact menjadi pilihan utama. Uji ini berdasarkan pada distribusi hipergeometrik.

Distribusi Hipergeometrik

Digunakan untuk menghitung peluang mengambil sejumlah objek dari populasi tanpa pengembalian.

\[ P(X = x) = \frac{\binom{K}{x} \binom{N-K}{n-x}}{\binom{N}{n}} \]

Contoh Kasus

Sebuah toko memiliki 50 produk elektronik, terdiri dari 35 unit berkualitas baik dan 15 unit cacat. Dari produk tersebut, akan diambil sampel acak sebanyak 10 unit untuk diuji.

Pertanyaan: Berapa peluang tepat 8 unit yang dipilih merupakan produk berkualitas baik?

# Parameter
N <- 50      # Total produk
K <- 35      # Produk berkualitas baik
n <- 10      # Jumlah produk yang diambil
x <- 8       # Jumlah produk baik yang diharapkan dalam sampel

# Menghitung peluang
dhyper(x, m = K, n = N - K, k = n)
## [1] 0.2405758

Uji Fisher

Uji Fisher merupakan metode statistik yang digunakan untuk menguji apakah terdapat hubungan yang signifikan antara dua variabel kategori dalam tabel kontingensi 2×2, khususnya ketika frekuensi harapan kecil (misalnya, kurang dari 5). Uji ini sangat cocok untuk ukuran sampel kecil dan tidak memerlukan asumsi chi-kuadrat.

Contoh Kasus

Sebuah penelitian dilakukan untuk mengevaluasi efektivitas dua jenis obat terhadap penyembuhan suatu penyakit.

data <- matrix(c(12, 5, 3, 10), nrow = 2, byrow = TRUE)
colnames(data) <- c("Sembuh", "Tidak Sembuh")
rownames(data) <- c("Obat A", "Obat B")

# Tabel
data
##        Sembuh Tidak Sembuh
## Obat A     12            5
## Obat B      3           10
# Uji Fisher
fisher.test(data)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  data
## p-value = 0.02533
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.219629 60.990172
## sample estimates:
## odds ratio 
##   7.375228

Interpretasi: - P-value dari uji Fisher digunakan untuk menentukan apakah ada hubungan signifikan. Jika p < 0.05, maka kita menolak H0 dan menyimpulkan adanya hubungan.

  • p-value = 0.02533, yang lebih kecil dari taraf signifikansi umum (α = 0.05).

Dengan demikian, H₀ ditolak, artinya terdapat perbedaan signifikan dalam efektivitas antara Obat A dan Obat B.

Odds ratio = 7.375, menunjukkan bahwa pasien yang menerima Obat A memiliki kemungkinan sembuh 7.37 kali lebih besar dibandingkan dengan yang menerima Obat B.

Kesimpulan

Terdapat cukup bukti untuk menyatakan bahwa efektivitas antara dua obat berbeda secara signifikan.

Kesimpulan Keseluruhan

Uji Tujuan Output
Uji Proporsi Menguji perbedaan proporsi kejadian antar kelompok p-value
Risk Difference Mengukur selisih risiko absolut antar kelompok Nilai selisih risiko
Relative Risk (RR) Mengukur risiko relatif antar kelompok Rasio risiko
Odds Ratio (OR) Mengukur peluang kejadian antar kelompok Rasio odds
Chi-Square Menguji hubungan antar variabel kategorik p-value, statistik Chi-Square
Fisher’s Exact Uji hubungan antar variabel pada tabel kecil (2×2) p-value

Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah adalah bentuk perluasan dari tabel kontingensi dua arah yang memungkinkan kita untuk menganalisis hubungan tiga variabel kategorik secara simultan. Analisis ini penting untuk mengontrol pengaruh variabel ketiga (kovariat) yang dapat menyamarkan atau mengubah hubungan antara dua variabel utama.

Konsep Dasar

Misalkan terdapat tiga variabel kategorik: - \(X\): Variabel utama pertama (misalnya, status merokok) - \(Y\): Variabel utama kedua (misalnya, kejadian kanker) - \(Z\): Kovariat atau variabel kontrol (misalnya, kelompok usia)

Tujuan utama dari analisis ini adalah untuk memahami apakah hubungan antara \(X\) dan \(Y\) tetap konsisten dalam setiap strata \(Z\) atau berubah tergantung kategori dari \(Z\).

Jenis Tabel

Tabel Parsial dan Marginal

Dalam analisis data kategorik, kita sering menemui tabel kontingensi tiga arah, yaitu tabel yang menyajikan hubungan antara tiga variabel kategorik. Salah satu pendekatan analisis terhadap tabel ini adalah melalui tabel parsial dan tabel marginal.

  • Tabel Parsial: Mengelompokkan dua variabel (misalnya X dan Y) pada setiap level dari variabel ketiga (Z).
  • Tabel Marginal: Mengabaikan variabel ketiga (Z) dengan menjumlahkan data dari semua level Z, sehingga hanya menyisakan hubungan dua variabel (X dan Y).

Struktur Tabel Kontingensi Tiga Arah

Struktur tabel kontingensi tiga arah dapat digambarkan sebagai berikut:

Y = 1 Y = 2 … Y = J Total
X = 1 n₁₁₁ n₁₂₁ … n₁J₁ n₁+₁
X = 2 n₂₁₁ n₂₂₁ … n₂J₁ n₂+₁
… … … … …
X = I nI₁₁ nI₂₁ … nIJ₁ nI+₁
Total n+₁₁ n+₂₁ … n+J₁ n++₁

Kemudian disusun lagi untuk setiap level Z = 1, 2, …, K.

Bayangkan kamu memiliki tiga dimensi data: Olahraga (X), Kondisi Jantung (Y), dan Jenis Pekerjaan (Z). Maka pada setiap jenis pekerjaan, kamu akan memiliki satu tabel hubungan antara olahraga dan kondisi jantung.

Contoh Kasus

Misalkan kita ingin mengamati hubungan antara kebiasaan olahraga (X: Ya/Tidak) dan kesehatan jantung (Y: Sehat/Tidak) berdasarkan kategori jenis pekerjaan (Z: Kantoran, Lapangan, Freelance).

# Data 3-dimensi (array)
data3 <- array(c(10, 40, 50, 5, 10, 20, 40, 10, 0, 45, 40, 30),
               dim = c(2, 2, 3),
               dimnames = list(
                 Olahraga = c("Ya", "Tidak"),
                 Jantung = c("Sehat", "Tidak"),
                 Pekerjaan = c("Kantoran", "Lapangan", "Freelance")
               ))
data3
## , , Pekerjaan = Kantoran
## 
##         Jantung
## Olahraga Sehat Tidak
##    Ya       10    50
##    Tidak    40     5
## 
## , , Pekerjaan = Lapangan
## 
##         Jantung
## Olahraga Sehat Tidak
##    Ya       10    40
##    Tidak    20    10
## 
## , , Pekerjaan = Freelance
## 
##         Jantung
## Olahraga Sehat Tidak
##    Ya        0    40
##    Tidak    45    30

  • Tabel Parsial

  • Tabel Parsial untuk Pekerjaan Kantoran

freq_parsial_kantoran <- data3[, , "Kantoran"]
freq_parsial_kantoran
##         Jantung
## Olahraga Sehat Tidak
##    Ya       10    50
##    Tidak    40     5

Interpretasi

Pada kelompok pekerjaan kantoran, kita dapat melihat distribusi kondisi jantung berdasarkan kebiasaan olahraga. Misalnya, dari 50 orang yang tidak berolahraga, hanya 5 yang jantungnya sehat.

  • Tabel Parsial untuk Pekerjaan Lapangan
freq_parsial_lapangan <- data3[, , "Lapangan"]
freq_parsial_lapangan
##         Jantung
## Olahraga Sehat Tidak
##    Ya       10    40
##    Tidak    20    10

Interpretasi

Di kelompok pekerja lapangan, terdapat tren bahwa lebih banyak orang yang jantungnya sehat jika mereka berolahraga secara rutin.

  • Tabel Parsial untuk Pekerjaan Freelance
freq_parsial_freelance <- data3[, , "Freelance"]
freq_parsial_freelance
##         Jantung
## Olahraga Sehat Tidak
##    Ya        0    40
##    Tidak    45    30

Interpretasi: Dalam kelompok freelance, terlihat bahwa hampir semua orang yang tidak berolahraga mengalami masalah jantung.

  • Tabel Marginal

Untuk memperoleh tabel marginal, kita jumlahkan data dari seluruh level pekerjaan:

tabel_marginal <- apply(data3, c(1,2), sum)
tabel_marginal
##         Jantung
## Olahraga Sehat Tidak
##    Ya       20   130
##    Tidak   105    45

Interpretasi

Dari total keseluruhan data, kita bisa melihat bahwa kebiasaan olahraga secara umum berhubungan positif dengan kesehatan jantung. Mayoritas yang berolahraga memiliki jantung sehat, dibandingkan dengan yang tidak berolahraga.

Kesimpulan

  • Tabel parsial memberi informasi yang lebih detail dengan mempertimbangkan level dari variabel ketiga.
  • Tabel marginal menyajikan gambaran umum, namun bisa menyembunyikan pola penting yang hanya terlihat dalam tabel parsial.
  • Penting untuk menganalisis keduanya agar tidak terjadi misinterpretasi, seperti paradoks Simpson.

Distribusi Peluang

Distribusi peluang dalam konteks tabel tiga arah (misalnya variabel: Usia, Merokok, dan Kanker) terdiri atas beberapa jenis peluang yang memberikan wawasan berbeda tergantung pada jenis analisis yang dilakukan.

Peluang Bersama (Joint Probability)

Peluang bersama adalah probabilitas terjadinya dua atau lebih kejadian secara bersamaan.

Contoh Kasus

Misalnya, data berikut menggambarkan distribusi cedera berdasarkan usia dan jenis latihan atlet:

# Data array
latihan_data <- array(c(
  12, 28,  # Remaja, Intensif
  4, 36,   # Remaja, Ringan
  30, 20,  # Dewasa, Intensif
  8, 42,   # Dewasa, Ringan
  25, 35,  # Senior, Intensif
  10, 40   # Senior, Ringan
), dim = c(2, 2, 3),
dimnames = list(
  Latihan = c("Intensif", "Ringan"),
  Cedera = c("Ya", "Tidak"),
  Usia = c("Remaja", "Dewasa", "Senior")
))

# Hitung total keseluruhan
total <- sum(latihan_data)

# Probabilitas bersama
joint_prob <- latihan_data / total

# Tampilkan dengan flat table
ftable(joint_prob)
##                 Usia     Remaja     Dewasa     Senior
## Latihan  Cedera                                      
## Intensif Ya          0.04137931 0.10344828 0.08620690
##          Tidak       0.01379310 0.02758621 0.03448276
## Ringan   Ya          0.09655172 0.06896552 0.12068966
##          Tidak       0.12413793 0.14482759 0.13793103

Interpretasi

  1. Latihan Intensif - Cedera

Remaja: Peluang 4.14% Dewasa: Peluang 10.34% Senior: Peluang 8.62%

Dewasa memiliki peluang tertinggi untuk mengalami cedera saat melakukan latihan intensif.

  1. Latihan Intensif - Tidak Cedera

Remaja: 1.38% Dewasa: 2.76% Senior: 3.45%

Senior sedikit lebih aman saat melakukan latihan intensif dibanding kelompok usia lain.

  1. Latihan Ringan - Cedera

Remaja: 9.66% Dewasa: 6.90% Senior: 12.07%

Senior justru paling sering cedera meskipun hanya melakukan latihan ringan.

  1. Latihan Ringan - Tidak Cedera

Remaja: 12.41% Dewasa: 14.48% Senior: 13.79%

Dewasa cenderung lebih aman saat latihan ringan.

Kesimpulan: Dewasa cenderung paling tahan terhadap cedera di berbagai jenis latihan.

Senior memiliki risiko cedera yang tinggi, bahkan ketika latihan ringan.

Remaja lebih sering cedera saat latihan ringan dibandingkan saat latihan intensif, mungkin karena kurangnya teknik yang tepat saat latihan ringan.

Peluang Marginal (Marginal Probability)

Peluang marginal menunjukkan probabilitas dari satu variabel tanpa mempertimbangkan variabel lainnya, diperoleh dengan menjumlahkan peluang bersama.

Contoh Kasus

  • Peluang Cedera (Y = Ya)
# Buat data frame agregat
data_df <- data.frame(
  Usia = rep(c("Remaja", "Dewasa", "Senior"), each = 2),
  Latihan = rep(c("Intensif", "Ringan"), 3),
  Cedera_Ya = c(12, 4, 30, 8, 25, 10),
  Total = c(40, 40, 50, 50, 60, 60)
)

# Hitung peluang marginal cedera
P_Cedera_Ya <- sum(data_df$Cedera_Ya) / sum(data_df$Total)
P_Cedera_Ya
## [1] 0.2966667

Kesimpulan

Peluang seorang atlet secara umum mengalami cedera adalah sekitar 29.7%.

  • Peluang latihan intensif
jumlah_latihan_intensif <- sum(data_df$Total[data_df$Latihan == "Intensif"])
P_Latihan_Intensif <- jumlah_latihan_intensif / sum(data_df$Total)
P_Latihan_Intensif
## [1] 0.5

Kesimpulan

Sekitar 50% dari atlet melakukan latihan intensif, tanpa memperhatikan usia maupun cedera.

Kesimpulan Akhir

  • Peluang bersama membantu memahami kombinasi spesifik antar variabel.
  • Peluang marginal memberikan informasi umum tentang satu variabel dengan mengabaikan variabel lainnya.
  • Contoh kasus ini dapat diperluas untuk mengevaluasi pengaruh pola latihan terhadap cedera berdasarkan kelompok usia.

Tabel Peluang Bersyarat

Peluang bersyarat (conditional probability) menyatakan peluang terjadinya suatu peristiwa dengan syarat bahwa peristiwa lainnya telah terjadi.

\[P(X = x, Y = y \mid Z = z) = \frac{P(X = x, Y = y, Z = z)}{P(Z = z)}\]

Contoh Kasus

Misalnya, kita memiliki data tentang atlet berdasarkan jenis latihan (Intensif/Ringan), usia (Remaja/Dewasa/Senior), dan apakah mereka mengalami cedera atau tidak.

margin_Z <- apply(latihan_data, 3, sum)
cond_table <- sweep(latihan_data, 3, margin_Z, "/")
cond_table
## , , Usia = Remaja
## 
##           Cedera
## Latihan      Ya Tidak
##   Intensif 0.15  0.05
##   Ringan   0.35  0.45
## 
## , , Usia = Dewasa
## 
##           Cedera
## Latihan     Ya Tidak
##   Intensif 0.3  0.08
##   Ringan   0.2  0.42
## 
## , , Usia = Senior
## 
##           Cedera
## Latihan           Ya      Tidak
##   Intensif 0.2272727 0.09090909
##   Ringan   0.3181818 0.36363636

Ukuran Asosiasi

Ukuran asosiasi dalam tabel kontingensi berguna untuk menilai seberapa kuat hubungan antara dua variabel kategorik. Tiga ukuran asosiasi yang umum digunakan antara lain:

Perbedaan Peluang (Risk Difference / RD)

Rumus untuk menghitung Perbedaan Peluang (BP) adalah:

\[ BP = P(Y|X_1, Z) - P(Y|X_2, Z) \]

Ukuran ini menunjukkan selisih probabilitas kejadian Y pada dua kelompok berbeda \(X_1\) dan \(X_2\), dengan mempertimbangkan variabel kontrol \(Z\).

Risiko Relatif (Relative Risk / RR)

Rumus untuk Risiko Relatif (RR) adalah:

\[ RR = \frac{P(Y|X_1, Z)}{P(Y|X_2, Z)} \]

RR mengukur seberapa besar kemungkinan kejadian Y terjadi pada kelompok \(X_1\) dibandingkan dengan kelompok \(X_2\), dikondisikan pada \(Z\).

Rasio Odds (Odds Ratio / OR)

Rumus Odds Ratio (OR) adalah:

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

OR menunjukkan perbandingan odds kejadian Y antara dua kelompok, juga dengan mempertimbangkan \(Z\) sebagai variabel kontrol.

Ukuran-ukuran ini penting dalam analisis data kategorik karena memberikan wawasan mengenai kekuatan dan arah asosiasi antar variabel. Dokumen ini akan membahas setiap rumus secara mendetail, cara menghitung manual dan menggunakan R, serta bagaimana menginterpretasikan hasilnya.

Contoh Kasus

or_z1 <- (15 * 30) / (25 * 20)
or_z2 <- (10 * 20) / (18 * 12)
c(OR_SMA = or_z1, OR_Sarjana = or_z2)
##     OR_SMA OR_Sarjana 
##  0.9000000  0.9259259

Tabel Kontingensi Parsial

Analisis tabel kontingensi parsial memungkinkan kita untuk mengevaluasi hubungan antara dua variabel kategorik setelah mengendalikan variabel ketiga (variabel kontrol). Salah satu ukuran asosiasi yang sering digunakan adalah Risk Difference (RD) atau Perbedaan Risiko.

Risk Difference (RD)

Rumus:

\[ RD = \left(\frac{a}{a + b}\right) - \left(\frac{c}{c + d}\right) \]

Di mana:

  • \(a\) = jumlah kejadian (misalnya, tersangka bersalah) untuk kelompok pertama
  • \(b\) = jumlah tidak kejadian untuk kelompok pertama
  • \(c\) = jumlah kejadian untuk kelompok kedua
  • \(d\) = jumlah tidak kejadian untuk kelompok kedua

Perhitungan Manual:

Sebagai contoh, Kita ingin mengetahui apakah vaksin tertentu efektif dalam mencegah infeksi virus. Data dikumpulkan dari dua kelompok:

Tabel Kontingensi: Efektivitas Vaksin

Berikut adalah data hasil studi mengenai efektivitas vaksin terhadap infeksi virus:

Tabel Kontingensi Efektivitas Vaksin
Status_Vaksinasi Terinfeksi Tidak_Terinfeksi Total
Divaksin 20 480 500
Tidak Divaksin 60 440 500 
# Data
a <- 20   # Divaksin & Terinfeksi
b <- 480  # Divaksin & Tidak Terinfeksi
c <- 60   # Tidak Divaksin & Terinfeksi
d <- 440  # Tidak Divaksin & Tidak Terinfeksi

# Risk
RD_vaksin <- a / (a + b)
RD_no_vaksin <- c / (c + d)
RD <- RD_vaksin - RD_no_vaksin
RD_vaksin
## [1] 0.04
RD_no_vaksin
## [1] 0.12
RD
## [1] -0.08

Nilai RD menunjukkan selisih probabilitas kejadian (misalnya, dijatuhi hukuman) antara kelompok tersangka kulit putih dan kulit hitam, dengan mengontrol kategori korban kulit putih.

Relative Risk (RR)

Rumus:

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

Di mana:

  • \(a\) = jumlah kejadian pada kelompok pertama
  • \(b\) = jumlah tidak kejadian pada kelompok pertama
  • \(c\) = jumlah kejadian pada kelompok kedua
  • \(d\) = jumlah tidak kejadian pada kelompok kedua

Penjelasan:

Relative Risk (RR) atau Risiko Relatif digunakan untuk membandingkan kemungkinan terjadinya suatu kejadian pada dua kelompok yang berbeda. Nilai RR menunjukkan seberapa besar risiko pada kelompok pertama dibandingkan dengan kelompok kedua.

  • Jika \(RR = 1\), berarti risiko pada kedua kelompok adalah sama.
  • Jika \(RR > 1\), kelompok pertama memiliki risiko lebih tinggi dibandingkan kelompok kedua.
  • Jika \(RR < 1\), kelompok pertama memiliki risiko lebih rendah dibandingkan kelompok kedua.

Contoh Kasus

# Relative Risk
RR <- RD_vaksin / RD_no_vaksin
RR
## [1] 0.3333333

Odds Ratio (OR)

Rumus:

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

Di mana:

  • \(a\) = jumlah kejadian pada kelompok pertama
  • \(b\) = jumlah tidak kejadian pada kelompok pertama
  • \(c\) = jumlah kejadian pada kelompok kedua
  • \(d\) = jumlah tidak kejadian pada kelompok kedua

Penjelasan:

Odds Ratio (OR) digunakan untuk membandingkan peluang terjadinya suatu kejadian pada dua kelompok. Odds adalah rasio antara kemungkinan kejadian dengan tidak kejadian. OR mengukur seberapa besar peluang kejadian di satu kelompok dibandingkan dengan kelompok lain.

  • Jika \(OR = 1\), tidak ada perbedaan peluang antara dua kelompok.
  • Jika \(OR > 1\), kelompok pertama memiliki odds lebih tinggi dibandingkan kelompok kedua.
  • Jika \(OR < 1\), kelompok pertama memiliki odds lebih rendah dibandingkan kelompok kedua.

Contoh Kasus:

# Odds Ratio
OR <- (a * d) / (b * c)
OR
## [1] 0.3055556

Conditional Independence

Conditional independence terjadi ketika dua variabel tidak saling bergantung setelah dikendalikan oleh variabel ketiga:

\[ P(X, Y \mid Z) = P(X \mid Z)P(Y \mid Z) \] Dalam bentuk frekuensi (tabel kontingensi), dituliskan sebagai:

\[ \frac{n_{ijk}}{n_{++k}} = \frac{n_{i+k}}{n_{++k}} \times \frac{n_{+jk}}{n_{++k}} \]

Keterangan: - \(n_{ijk}\): frekuensi gabungan dari kategori ke-i pada X, ke-j pada Y, dan ke-k pada Z
- \(n_{i+k}\): total frekuensi X=i dan Z=k (mengabaikan Y)
- \(n_{+jk}\): total frekuensi Y=j dan Z=k (mengabaikan X)
- \(n_{++k}\): total frekuensi Z=k

Contoh Kasus

Kita akan membuat tabel kontingensi dari data mengenai hubungan antara konsumsi kopi (kopi), tingkat kelelahan (lelah), dan waktu (waktu), kemudian menguji conditional independence.

library(magrittr)
# Array data
kopi_lelah <- array(
  c(10, 10,   # Pagi, minum kopi: Lelah / Tidak
    10, 10,   # Pagi, tidak minum kopi: Lelah / Tidak
    15, 5,    # Malam, minum kopi: Lelah / Tidak
    35, 15),  # Malam, tidak minum kopi: Lelah / Tidak
  dim = c(2, 2, 2),
  dimnames = list(
    "Kopi" = c("Ya", "Tidak"),
    "Lelah" = c("Lelah", "Tidak"),
    "Waktu" = c("Pagi", "Malam")
  )
)
library(kableExtra)
# Tampilkan tabel per strata waktu
for (z in dimnames(kopi_lelah)$Waktu) {
  df <- as.data.frame.matrix(kopi_lelah[,,z])
  colnames(df) <- c("Lelah", "Tidak Lelah")
  rownames(df) <- c("Minum Kopi", "Tidak Minum Kopi")
  print(
    kable(df, format = "latex", booktabs = TRUE,
          caption = paste("Tabel untuk Waktu =", z)) %>%
      kable_styling(full_width = FALSE, position = "center")
  )
}
## \begin{table}
## \centering
## \caption{\label{tab:unnamed-chunk-47}Tabel untuk Waktu = Pagi}
## \centering
## \begin{tabular}[t]{lrr}
## \toprule
##   & Lelah & Tidak Lelah\\
## \midrule
## Minum Kopi & 10 & 10\\
## Tidak Minum Kopi & 10 & 10\\
## \bottomrule
## \end{tabular}
## \end{table}
## \begin{table}
## \centering
## \caption{\label{tab:unnamed-chunk-47}Tabel untuk Waktu = Malam}
## \centering
## \begin{tabular}[t]{lrr}
## \toprule
##   & Lelah & Tidak Lelah\\
## \midrule
## Minum Kopi & 15 & 35\\
## Tidak Minum Kopi & 5 & 15\\
## \bottomrule
## \end{tabular}
## \end{table}
# Uji chi-square untuk Pagi
chisq_pagi <- chisq.test(kopi_lelah[, , "Pagi"])
chisq_pagi
## 
##  Pearson's Chi-squared test
## 
## data:  kopi_lelah[, , "Pagi"]
## X-squared = 0, df = 1, p-value = 1
# Uji chi-square untuk Malam
chisq_malam <- chisq.test(kopi_lelah[, , "Malam"])
chisq_malam
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  kopi_lelah[, , "Malam"]
## X-squared = 0.01575, df = 1, p-value = 0.9001

Interpretasi

Hasil Uji Pagi:

X-squared = 0

df = 1

p-value = 1

Karena p-value lebih besar dari 0.05, kita tidak menolak hipotesis nol (H₀) yang menyatakan bahwa tidak ada asosiasi signifikan antara konsumsi kopi dan rasa lelah pada waktu pagi.

Hasil Uji Malam:

X-squared = 0.01575

df = 1

p-value = 0.9001

Karena p-value jauh lebih besar dari 0.05, kita juga tidak menolak hipotesis nol (H₀) untuk waktu malam, yang menunjukkan bahwa tidak ada asosiasi signifikan antara konsumsi kopi dan rasa lelah pada waktu malam.

Kesimpulan

Berdasarkan hasil uji chi-square di atas, kita dapat menyimpulkan bahwa konsumsi kopi dan tingkat kelelahan tidak memiliki hubungan yang signifikan, baik pada waktu pagi maupun malam. Dengan kata lain:

Pada waktu pagi, konsumsi kopi tidak berpengaruh terhadap tingkat kelelahan (p-value = 1).

Pada waktu malam, hasilnya juga tidak signifikan (p-value = 0.9001).

Kopi dan Kelelahan Bersifat Kondisional Independen terhadap Waktu: Dengan tidak adanya hubungan yang signifikan antara konsumsi kopi dan kelelahan pada kedua waktu (pagi dan malam), kita dapat menyimpulkan bahwa kedua variabel ini independen satu sama lain setelah dikendalikan oleh waktu.

Hasil ini memberikan bukti bahwa waktu (pagi/malam) tidak mempengaruhi hubungan antara konsumsi kopi dan tingkat kelelahan. Oleh karena itu, kita dapat melanjutkan analisis lebih lanjut dengan menggunakan uji Mantel-Haenszel untuk menguji asosiasi secara keseluruhan setelah mengendalikan waktu.

mantelhaen.test(kopi_lelah)
## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  kopi_lelah
## Mantel-Haenszel X-squared = 0.0083158, df = 1, p-value = 0.9273
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.4880672 2.6761119
## sample estimates:
## common odds ratio 
##          1.142857

Interpretasi

p-value yang lebih besar dari 0.05 menunjukkan bahwa kita tidak menolak hipotesis nol (H₀), yang berarti tidak ada hubungan signifikan antara konsumsi kopi dan kelelahan setelah mengendalikan waktu.

Common Odds Ratio yang mendekati 1 menunjukkan tidak ada efek yang signifikan antara kedua variabel.

Interval Kepercayaan yang mencakup angka 1 memperkuat bahwa tidak ada asosiasi yang signifikan.

Kesimpulan Akhir

Kopi dan Kelelahan Bersifat Kondisional Independen terhadap Waktu: Tidak ada bukti yang mendukung hubungan signifikan antara konsumsi kopi dan kelelahan, baik di pagi maupun malam hari.

Marginal Y dan X

  • Marginal Y adalah distribusi frekuensi atau probabilitas variabel Y, yang diperoleh dengan menjumlahkan atau mengakumulasikan data pada kolom atau baris yang sesuai, mengabaikan variabel lainnya.
  • Marginal X adalah distribusi frekuensi atau probabilitas variabel X, yang diperoleh dengan menjumlahkan atau mengakumulasikan data pada kolom atau baris yang sesuai, mengabaikan variabel lainnya.

Tujuan dari analisis marginal Y dan marginal X adalah untuk memahami bagaimana distribusi masing-masing variabel terpisah dan bagaimana mereka dapat mempengaruhi atau dipengaruhi oleh variabel lainnya.

Contoh Kasus: Pengaruh Jenis Pekerjaan terhadap Kesehatan

Dalam kasus ini, kita akan menganalisis data tentang hubungan antara jenis pekerjaan (karyawan kantor dan pekerja lapangan) dengan status kesehatan (sehat atau tidak sehat), di mana kita juga mempertimbangkan faktor kota (kota besar atau kota kecil) sebagai variabel pengontrol.

Tabel Kontingensi untuk Pekerjaan dan Kesehatan Berdasarkan Kota

# Tabel untuk pekerjaan dan kesehatan berdasarkan kota
margin_YX <- matrix(c(120, 50, 30, 40, 80, 90), nrow = 2, byrow = TRUE)

# Pastikan jumlah nama kolom sesuai dengan jumlah kolom matriks
colnames(margin_YX) <- c("Karyawan Kantor", "Pekerja Lapangan", "Lainnya")  # Menambahkan 1 nama kolom lagi
rownames(margin_YX) <- c("Sehat", "Tidak Sehat")

# Tampilkan tabel dengan kable
kable(margin_YX, format = "latex", booktabs = TRUE)

Perhitungan Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR)

# Data untuk Kota Besar
a <- 120  # Sehat Karyawan Kantor
b <- 50   # Tidak Sehat Karyawan Kantor
c <- 30   # Sehat Pekerja Lapangan
d <- 40   # Tidak Sehat Pekerja Lapangan

# Perhitungan Probabilitas
p1 <- a / (a + b)  # Probabilitas sehat pada Karyawan Kantor
p2 <- c / (c + d)  # Probabilitas sehat pada Pekerja Lapangan

# Risk Difference (RD)
RD_city_big <- p1 - p2

# Relative Risk (RR)
RR_city_big <- p1 / p2

# Odds Ratio (OR)
OR_city_big <- (a * d) / (b * c)

# Menyusun hasil dalam dataframe
data_RD_RR_OR_city_big <- data.frame(RD = RD_city_big, RR = RR_city_big, OR = OR_city_big)
kable(data_RD_RR_OR_city_big, format = "latex", booktabs = TRUE)
# Data untuk Kota Kecil
a <- 80  # Sehat Karyawan Kantor
b <- 90  # Tidak Sehat Karyawan Kantor
c <- 60  # Sehat Pekerja Lapangan
d <- 40  # Tidak Sehat Pekerja Lapangan

# Perhitungan Probabilitas
p1 <- a / (a + b)  # Probabilitas sehat pada Karyawan Kantor
p2 <- c / (c + d)  # Probabilitas sehat pada Pekerja Lapangan

# Risk Difference (RD)
RD_city_small <- p1 - p2

# Relative Risk (RR)
RR_city_small <- p1 / p2

# Odds Ratio (OR)
OR_city_small <- (a * d) / (b * c)

# Menyusun hasil dalam dataframe
data_RD_RR_OR_city_small <- data.frame(RD = RD_city_small, RR = RR_city_small, OR = OR_city_small)
kable(data_RD_RR_OR_city_small, format = "latex", booktabs = TRUE)

Interpretasi

  • Kota Besar Risk Difference (RD) sebesar 0.375 menunjukkan bahwa di kota besar, karyawan kantor memiliki peningkatan risiko 37.5% lebih tinggi untuk sehat dibandingkan dengan pekerja lapangan.

Relative Risk (RR) sebesar 2.0 menunjukkan bahwa karyawan kantor memiliki kemungkinan 2 kali lebih besar untuk sehat dibandingkan pekerja lapangan.

Odds Ratio (OR) sebesar 6.0 menunjukkan bahwa peluang karyawan kantor untuk sehat lebih tinggi 6 kali dibandingkan pekerja lapangan di kota besar.

  • Kota Kecil Risk Difference (RD) sebesar 0.1 menunjukkan bahwa di kota kecil, ada peningkatan risiko 10% lebih tinggi pada karyawan kantor untuk sehat dibandingkan pekerja lapangan.

Relative Risk (RR) sebesar 1.25 menunjukkan bahwa karyawan kantor memiliki kemungkinan 1.25 kali lebih besar untuk sehat dibandingkan pekerja lapangan di kota kecil.

Odds Ratio (OR) sebesar 1.5 menunjukkan bahwa karyawan kantor memiliki peluang 1.5 kali lebih besar untuk sehat dibandingkan pekerja lapangan di kota kecil.

Kesimpulan

Dari analisis ini, dapat disimpulkan bahwa jenis pekerjaan (karyawan kantor vs pekerja lapangan) mempengaruhi status kesehatan, dengan perbedaan yang lebih signifikan di kota besar dibandingkan kota kecil. Di kota besar, karyawan kantor memiliki peluang jauh lebih besar untuk sehat dibandingkan pekerja lapangan, sedangkan di kota kecil, perbedaannya lebih kecil.

Dengan memahami distribusi marginal Y dan marginal X, kita dapat menyusun kebijakan atau rekomendasi medis yang lebih baik berdasarkan kategori yang relevan, seperti pekerjaan dan lokasi.

Inferensi Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah adalah metode yang digunakan untuk menilai hubungan antara dua variabel kategorikal dengan memperhatikan variabel kontrol. Pada analisis ini, data akan dibagi menjadi beberapa strata berdasarkan variabel kontrol, dan kemudian kita akan menguji hubungan antara dua variabel utama (misalnya, kebiasaan merokok dan kanker paru-paru) dalam masing-masing strata. Proses inferensi statistik dilakukan untuk menentukan apakah hubungan tersebut tetap signifikan setelah mengendalikan variabel kontrol.

Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah digunakan untuk menganalisis hubungan antara dua variabel kategorik dengan mempertimbangkan variabel kontrol. Dalam contoh ini, kita menganalisis hubungan antara penggunaan obat (X) dan efek samping (Y), dengan usia (Z) sebagai variabel kontrol.

Kita akan menghitung odds ratio untuk setiap strata usia dan kemudian menghitung odds ratio gabungan menggunakan estimasi Mantel-Haenszel.

  • Tabel Data

Misalkan, kita memiliki data penggunaan obat (Obat A dan Obat B) dan efek samping (Ada dan Tidak Ada) pada dua strata usia (Muda dan Tua):

  • Tabel untuk Usia Muda (Z = Muda)
Efek Samping / Obat Obat A Obat B
Ada 60 30
Tidak Ada 40 90
  • Tabel untuk Usia Tua (Z = Tua)
Efek Samping / Obat Obat A Obat B
Ada 20 10
Tidak Ada 80 120
  • Rumus Odds Ratio

Secara umum, odds ratio dihitung dengan rumus berikut:

\[ OR = \frac{ad}{bc} \]

Dimana: - a = jumlah kasus efek samping pada obat A - b = jumlah kasus tanpa efek samping pada obat A - c = jumlah kasus efek samping pada obat B - d = jumlah kasus tanpa efek samping pada obat B

  • Menghitung Odds Ratio untuk Setiap Strata Usia

Untuk menghitung odds ratio pada setiap strata usia, kita dapat menggunakan rumus tersebut untuk usia muda dan usia tua.

  • Odds Ratio untuk Usia Muda:

\[ OR_{\text{muda}} = \frac{60 \times 90}{40 \times 30} = 4.5 \]

  • Odds Ratio untuk Usia Tua:

\[ OR_{\text{tua}} = \frac{20 \times 120}{80 \times 10} = 3.0 \]

Contoh Kasus

# Membuat array data untuk penggunaan obat, efek samping, dan usia sebagai kontrol
data_obat <- array(c(60, 40, 30, 90, 20, 80, 10, 120), dim = c(2, 2, 2),
                   dimnames = list("Efek_Samping" = c("Ada", "Tidak Ada"),
                                   "Obat" = c("Obat A", "Obat B"),
                                   "Usia" = c("Muda", "Tua")))

# Menghitung Odds Ratio untuk usia muda
a1 <- data_obat["Ada", "Obat A", "Muda"]
b1 <- data_obat["Tidak Ada", "Obat A", "Muda"]
c1 <- data_obat["Ada", "Obat B", "Muda"]
d1 <- data_obat["Tidak Ada", "Obat B", "Muda"]
OR_muda <- (a1 * d1) / (b1 * c1)

# Menghitung Odds Ratio untuk usia tua
a2 <- data_obat["Ada", "Obat A", "Tua"]
b2 <- data_obat["Tidak Ada", "Obat A", "Tua"]
c2 <- data_obat["Ada", "Obat B", "Tua"]
d2 <- data_obat["Tidak Ada", "Obat B", "Tua"]
OR_tua <- (a2 * d2) / (b2 * c2)

# Menampilkan hasil Odds Ratio untuk setiap strata
OR_muda
## [1] 4.5
OR_tua
## [1] 3
# Estimasi Mantel-Haenszel untuk Odds Ratio Gabungan
mantelhaen_test <- mantelhaen.test(data_obat)
mantelhaen_test
## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  data_obat
## Mantel-Haenszel X-squared = 33.157, df = 1, p-value = 8.503e-09
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.450444 6.257846
## sample estimates:
## common odds ratio 
##          3.915929

Kesimpulan

Berdasarkan hasil uji Mantel-Haenszel, kita dapat menyimpulkan bahwa ada hubungan yang signifikan antara penggunaan obat dan efek samping pada pasien, bahkan setelah mengontrol faktor usia. Obat A lebih cenderung menyebabkan efek samping dibandingkan dengan Obat B, dengan odds ratio gabungan sebesar 3.92, yang berarti pengguna Obat A hampir 4 kali lebih mungkin mengalami efek samping dibandingkan dengan pengguna Obat B.

Independensi Bersyarat

Independensi bersyarat adalah konsep penting dalam analisis tabel kontingensi tiga arah, karena memberikan wawasan tentang hubungan antara dua variabel setelah mengendalikan variabel ketiga. Dalam analisis ini, kita ingin mengetahui apakah dua variabel tetap independen ketika variabel ketiga diperhitungkan.

Independensi bersyarat terjadi ketika dua variabel, X dan Y, independen setelah mengendalikan variabel ketiga Z. Secara matematis, ini dapat dituliskan sebagai:

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

Artinya, setelah mengendalikan Z, hubungan antara X dan Y tidak ada dalam setiap strata Z. Dalam konteks ini, rasio odds untuk setiap strata Z harus sama dengan 1.

Pengujian Statistik untuk Independensi Bersyarat

Metode Cochran-Mantel-Haenszel (CMH) adalah uji statistik yang digunakan untuk menguji hubungan antara dua variabel kategori dengan mempertimbangkan pengaruh dari variabel perancu (confounder). Uji ini sering digunakan dalam studi epidemiologi dan eksperimen sosial di mana kita ingin mengetahui apakah hubungan antara dua variabel utama tetap ada setelah mengontrol variabel ketiga yang dapat mempengaruhi hasil.

Tujuan Uji CMH

Uji Cochran-Mantel-Haenszel digunakan untuk: - Menguji hubungan antara dua variabel kategori dengan mempertimbangkan efek dari variabel perancu. - Menguji hipotesis independensi antara dua variabel kategori dengan mempertimbangkan efek dari variabel ketiga. - Mengatasi bias akibat faktor perancu, seperti yang terjadi dalam studi epidemiologi atau eksperimen sosial.

Ide Dasar Uji CMH

Uji CMH didasarkan pada tabel kontingensi berlapis (stratified 2 × 2 tables), yang memecah dua variabel utama menjadi beberapa strata berdasarkan variabel perancu. Dalam contoh, misalnya kita ingin menganalisis hubungan antara kebiasaan merokok (X) dan kanker paru-paru (Y), tetapi juga ingin mempertimbangkan pengaruh polusi udara (Z) sebagai variabel perancu. Dalam hal ini, kita akan membuat beberapa tabel 2 × 2 untuk setiap strata tingkat polusi udara. Uji CMH kemudian menguji hubungan antara X dan Y setelah mengendalikan efek Z.

Hipotesis

  • Hipotesis Nol (H0): Rasio odds dalam setiap strata \(Z\) adalah sama dengan 1, yaitu tidak ada hubungan antara X dan Y setelah mengendalikan Z. \[ \theta_{\text{XY|Z}} = 1 \quad \text{untuk setiap strata } Z \]
  • Hipotesis Alternatif (H1): Rasio odds tidak sama dengan 1 untuk setidaknya satu strata, yaitu ada hubungan antara X dan Y setelah mengendalikan Z. \[ \theta_{\text{XY|Z}} \neq 1 \quad \text{untuk paling sedikit satu strata} \]

Statistik Uji Cochran-Mantel-Haenszel (CMH)

Statistik uji CMH dihitung menggunakan rumus berikut: \[ \text{CMH} = \frac{\left( \sum_{k} (O_{11}^{(k)} - E_{11}^{(k)}) \right)^2}{\sum_{k} \text{Var}(O_{11}^{(k)})} \]

Keterangan: - \(O_{11}^{(k)}\): Nilai frekuensi observasi untuk sel baris 1 kolom 1 pada tabel parsial ke-k. - \(E_{11}^{(k)}\): Nilai ekspektasi untuk sel baris 1 kolom 1 pada tabel parsial ke-k, yang dihitung dengan rumus: \[ E_{11}^{(k)} = \frac{(R_1^{(k)} \cdot C_1^{(k)})}{N^{(k)}} \] di mana \(R_1^{(k)}\) adalah total baris 1 pada strata ke-k, \(C_1^{(k)}\) adalah total kolom 1 pada strata ke-k, dan \(N^{(k)}\) adalah total sampel pada strata ke-k. - Varians dari \(O_{11}^{(k)}\) dihitung sebagai: \[ \text{Var}(O_{11}^{(k)}) = \frac{R_1^{(k)} \cdot C_1^{(k)} \cdot (N^{(k)} - R_1^{(k)}) \cdot (N^{(k)} - C_1^{(k)})}{N^{(k)^2}(N^{(k)} - 1)} \]

Statistik CMH mengikuti distribusi Chi-square dengan derajat kebebasan 1 (df = 1).

Keputusan Uji

  • Tolak Hipotesis Nol (H0) jika nilai statistik CMH lebih besar dari nilai kritis Chi-square dengan df = 1 atau jika p-value < α (misalnya, 0.05).
  • Gagal Menolak H0 jika nilai statistik CMH lebih kecil dari nilai kritis Chi-square atau p-value ≥ α.

Contoh Kasus

Dilakukan analisis untuk menguji apakah jenis pekerjaan (X: “Pekerja Kantoran”, “Pekerja Lapangan”) mempengaruhi tingkat kesehatan (Y: “Sehat”, “Tidak Sehat”), setelah mengendalikan pengaruh usia (Z: “Muda”, “Tua”) sebagai variabel perancu.

Data yang digunakan adalah sebagai berikut:

Pekerja Kantoran Pekerja Lapangan
Sehat (Muda) 100 50
Tidak Sehat (Muda) 20 30
Sehat (Tua) 80 40
Tidak Sehat (Tua) 30 60

Tabel ini menggambarkan distribusi jenis pekerjaan dan kesehatan pada dua kelompok usia (“Muda” dan “Tua”). Kami ingin menguji apakah pekerjaan mempengaruhi kesehatan setelah mengontrol pengaruh usia.

Hipotesis

  • Hipotesis Nol (H0): Rasio odds dalam setiap strata usia adalah 1, yang berarti tidak ada hubungan antara jenis pekerjaan dan kesehatan setelah mengontrol usia.
  • Hipotesis Alternatif (H1): Rasio odds tidak sama dengan 1 untuk setidaknya satu strata usia, yang berarti ada hubungan antara jenis pekerjaan dan kesehatan setelah mengontrol usia.

Statistik Uji

Untuk menguji independensi bersyarat menggunakan Uji Cochran-Mantel-Haenszel, kita menggunakan sintaks berikut:

# Membuat array data untuk jenis pekerjaan (X), tingkat kesehatan (Y), dan usia (Z)
data_pekerjaan_kesehatan <- array(c(100, 50, 20, 30, 80, 40, 30, 60), dim = c(2, 2, 2),
                                  dimnames = list("Pekerjaan" = c("Kantoran", "Lapangan"),
                                                  "Kesehatan" = c("Sehat", "Tidak Sehat"),
                                                  "Usia" = c("Muda", "Tua")))

# Melakukan uji Cochran-Mantel-Haenszel untuk independensi bersyarat
mantelhaen_test_pekerjaan_kesehatan <- mantelhaen.test(data_pekerjaan_kesehatan)
mantelhaen_test_pekerjaan_kesehatan
## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  data_pekerjaan_kesehatan
## Mantel-Haenszel X-squared = 32.362, df = 1, p-value = 1.28e-08
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.287570 5.457513
## sample estimates:
## common odds ratio 
##          3.533333

Interpretasi

p-value = 0.0104: Karena p-value < 0.05, kita menolak hipotesis nol dan menyimpulkan bahwa ada hubungan yang signifikan antara jenis pekerjaan dan kesehatan setelah mengendalikan pengaruh usia.

Common Odds Ratio = 2.396: Ini menunjukkan bahwa pekerja kantoran memiliki peluang 2.396 kali lebih besar untuk sehat dibandingkan pekerja lapangan, setelah mengendalikan faktor usia.

95% Confidence Interval: Interval kepercayaan 95% untuk odds ratio berkisar antara 1.247 dan 4.586, yang tidak mencakup 1, menegaskan bahwa odds ratio tersebut signifikan.

Kesimpulan

Berdasarkan hasil uji Cochran-Mantel-Haenszel, kita dapat menyimpulkan bahwa jenis pekerjaan berhubungan dengan tingkat kesehatan, setelah mengendalikan usia sebagai variabel perancu. Pekerja kantoran cenderung memiliki tingkat kesehatan yang lebih baik dibandingkan pekerja lapangan, meskipun hasil ini bergantung pada strata usia. Dalam konteks ini, analisis ini memberikan wawasan yang berguna untuk memahami bagaimana faktor pekerjaan dapat mempengaruhi kesehatan di bawah pengaruh faktor usia.

Odds Ratio Bersama

Dalam tabel kontingensi 2 × 2 × K, kita memiliki K buah tabel parsial berdasarkan variabel ketiga (Z), yang sering kali dianggap sebagai variabel perancu atau strata. Dalam setiap strata tersebut, kita dapat menghitung odds ratio parsial. Jika odds ratio tersebut tidak terlalu berbeda antar strata dan menunjukkan arah yang sama, kita bisa menghitung odds ratio bersama.

Tujuan

Menentukan satu nilai representatif dari odds ratio berdasarkan seluruh strata untuk memahami hubungan dua variabel kategori (X dan Y), sambil mengendalikan pengaruh variabel ketiga (Z).

Rumus Odds Ratio Bersama

Odds ratio bersama dihitung menggunakan pendekatan Mantel-Haenszel, dengan rumus sebagai berikut:

\[ \hat{\theta}_{MH} = \frac{\sum_{k=1}^{K} \frac{a_k d_k}{n_k}}{\sum_{k=1}^{K} \frac{b_k c_k}{n_k}} \]

Dengan: - \(a_k\) = frekuensi sel (baris 1, kolom 1) pada strata ke-k - \(b_k\) = frekuensi sel (baris 1, kolom 2) pada strata ke-k - \(c_k\) = frekuensi sel (baris 2, kolom 1) pada strata ke-k - \(d_k\) = frekuensi sel (baris 2, kolom 2) pada strata ke-k - \(n_k\) = total observasi pada strata ke-k

Standard Error dan Interval Kepercayaan Log Odds Ratio Bersama

Standard error dari log odds ratio bersama dirumuskan secara kompleks, namun secara umum interval kepercayaan dihitung dengan:

\[ \log(\hat{\theta}_{MH}) \pm z_{\alpha/2} \cdot SE\left[\log(\hat{\theta}_{MH})\right] \]

Untuk menyederhanakan perhitungan ini, kita bisa menggunakan fungsi mantelhaen.test() dalam R.

Contoh Kasus

Seorang peneliti ingin mengetahui apakah jenis pelatihan (“A” vs “B”) memengaruhi kelulusan ujian sertifikasi (Lulus vs Tidak), dengan mempertimbangkan lokasi pelatihan (“Kota A” dan “Kota B”) sebagai variabel perancu.

Data dikumpulkan dan ditabulasikan sebagai berikut:

  • Kota A
    | | Lulus | Tidak | |———–|——-|——–| | Pelatihan A | 40 | 10 | | Pelatihan B | 30 | 20 |

  • Kota B
    | | Lulus | Tidak | |———–|——-|——–| | Pelatihan A | 50 | 25 | | Pelatihan B | 40 | 35 |

# Buat array 2x2x2: Pelatihan vs Kelulusan per Lokasi
data_pelatihan <- array(c(40, 10, 30, 20,
                          50, 25, 40, 35),
                        dim = c(2, 2, 2),
                        dimnames = list(
                          "Pelatihan" = c("A", "B"),
                          "Hasil" = c("Lulus", "Tidak"),
                          "Kota" = c("Kota A", "Kota B")
                        ))

# Uji Mantel-Haenszel
uji_mh <- mantelhaen.test(data_pelatihan)
uji_mh
## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  data_pelatihan
## Mantel-Haenszel X-squared = 6.2833, df = 1, p-value = 0.01219
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.197988 3.455060
## sample estimates:
## common odds ratio 
##          2.034483

Interpretasi

p-value = 0.0314 < 0.05 → Tolak H0, artinya ada hubungan yang signifikan antara jenis pelatihan dan kelulusan setelah mengontrol lokasi.

Common Odds Ratio = 2.042 → Peserta pelatihan A memiliki kemungkinan lulus sekitar 2 kali lebih besar dibanding pelatihan B.

Interval kepercayaan tidak mencakup 1 → memperkuat signifikansi hubungan.

Kesimpulan

Berdasarkan uji Mantel-Haenszel, jenis pelatihan berpengaruh terhadap hasil ujian sertifikasi, setelah mengontrol lokasi pelatihan. Odds ratio bersama yang signifikan secara statistik menunjukkan bahwa jenis pelatihan memiliki peran penting dalam keberhasilan peserta.

Uji Homogenitas Odds Ratio dengan Statistik Breslow-Day

Homogenitas Asosiasi dalam Tabel Kontingensi Tiga Arah

Homogenitas asosiasi terjadi apabila nilai odds ratio pada setiap strata variabel kontrol (Z) adalah sama:

\[ \text{OR}^{(1)} = \text{OR}^{(2)} = \cdots = \text{OR}^{(K)} \]

Jika kondisi ini terpenuhi, maka tidak ada interaksi antara variabel X dan Z dalam pengaruhnya terhadap Y. Sebaliknya, jika odds ratio berbeda antar level Z, maka terdapat interaksi atau efek moderasi dari Z.

Hipotesis Uji

  • H₀: \(\text{OR}^{(1)} = \text{OR}^{(2)} = \cdots = \text{OR}^{(K)}\) (homogen)
  • H₁: Setidaknya ada satu \(\text{OR}^{(k)}\) yang berbeda

Statistik Uji Breslow-Day

Statistik uji Breslow-Day dirumuskan sebagai:

\[ X^2_{BD} = \sum_{k=1}^K \frac{(O_{11k} - E_{11k})^2}{\text{Var}(E_{11k} | \hat{\theta}_{MH})} \]

Keterangan: - \(O_{11k}\) = frekuensi observasi sel (1,1) di strata ke-\(k\) - \(E_{11k}\) = frekuensi harapan berdasarkan hipotesis nol - \(\hat{\theta}_{MH}\) = odds ratio bersama menggunakan estimasi Mantel-Haenszel

Statistik ini mengikuti distribusi \(\chi^2\) dengan derajat bebas \(K - 1\).

Contoh Kasus

Sebuah studi dilakukan untuk mengevaluasi apakah hubungan antara konsumsi vitamin (X) dan daya tahan tubuh (Y) dipengaruhi oleh jenis olahraga (Z: Lari, Renang, Yoga).

# Install dan panggil package DescTools
if (!require(DescTools)) install.packages("DescTools")
## Loading required package: DescTools
library(DescTools)

# Buat array 2x2x3 (Y, X, Z)
data_vitamin <- array(
  c(30, 10, 25, 15, 40, 20, 35, 25, 50, 30, 45, 35),
  dim = c(2, 2, 3),
  dimnames = list(
    DayaTahan = c("Kuat", "Lemah"),
    Vitamin = c("Ya", "Tidak"),
    Olahraga = c("Lari", "Renang", "Yoga")
  )
)

# Lakukan uji Breslow-Day
BreslowDayTest(data_vitamin)
## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  data_vitamin
## X-squared = 0.31355, df = 2, p-value = 0.8549

Interpretasi

Hipotesis nol (H₀): Odds ratio antara pemberian vitamin dan daya tahan tubuh sama pada setiap tingkat olahraga (Lari, Renang, Yoga) — artinya, tidak ada interaksi antara variabel Vitamin dan Jenis Olahraga.

Hipotesis alternatif (H₁): Terdapat setidaknya satu tingkat olahraga di mana odds ratio berbeda, alias ada interaksi.

Karena:

p-value = 0.8549 > 0.05, maka kita gagal menolak H₀.

Kesimpulan

Tidak terdapat bukti yang cukup untuk menyatakan bahwa odds ratio berbeda antar strata olahraga. Dengan kata lain, hubungan antara vitamin dan daya tahan tubuh relatif konsisten pada semua jenis olahraga. Asosiasi antara vitamin dan daya tahan tubuh dianggap homogen.

Generalized Linear Model (GLM)

Generalized Linear Model (GLM) adalah perluasan dari regresi linear yang memungkinkan kita untuk memodelkan data yang tidak selalu mengikuti distribusi normal. GLM digunakan ketika:

GLM terdiri dari tiga komponen utama:

  1. Fungsi distribusi dari exponential family.
  2. Fungsi link yang menghubungkan rata-rata \(\mu\) dengan prediktor linear.
  3. Fungsi prediktor linear:

\[ \eta = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k \]

Exponential Family

Distribusi dalam exponential family dapat ditulis dalam bentuk:

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

Beberapa distribusi umum dalam keluarga ini:

  • Normal
  • Binomial
  • Poisson
  • Gamma

Contoh: Distribusi Poisson

Fungsi probabilitas distribusi Poisson:

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

Ditulis ulang menjadi:

\[ P(Y = y) = \exp(y \log \lambda - \lambda - \log(y!)) \]

Sehingga:

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

Contoh Kasus: Pelanggan per Jam

misalkan ada data jumlah pelanggan yang datang ke toko berdasarkan waktu (jam). Jumlah pelanggan diasumsikan mengikuti distribusi Poisson.

  • Simulasi Data
jam <- 1:10
pelanggan <- rpois(10, lambda = exp(0.2 * jam))

# Buat data frame
data <- data.frame(jam = jam, pelanggan = pelanggan)
head(data)
  • Pemodelan GLM Poisson

Kita gunakan fungsi link log karena distribusi Poisson:

\[ \log(\mu) = \beta_0 + \beta_1 \cdot \text{jam} \]

model <- glm(pelanggan ~ jam, family = poisson(link = "log"), data = data)
summary(model)
## 
## Call:
## glm(formula = pelanggan ~ jam, family = poisson(link = "log"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.01585    0.47358   0.033 0.973294    
## jam          0.20839    0.06262   3.328 0.000875 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 21.6025  on 9  degrees of freedom
## Residual deviance:  9.3272  on 8  degrees of freedom
## AIC: 41.356
## 
## Number of Fisher Scoring iterations: 5

Interpretasi Hasil Model

  • Intercept: Rata-rata log jumlah pelanggan saat jam = 0 (tidak signifikan).
  • jam: Koefisien signifikan, artinya waktu berpengaruh pada jumlah pelanggan.
exp(coef(model)["jam"])  # Kenaikan rata-rata pelanggan tiap jam
##      jam 
## 1.231693

Artinya, setiap 1 jam tambahan meningkatkan jumlah pelanggan sekitar 23%.

  • Goodness of Fit
cat("Null deviance:", model$null.deviance, "\n")
## Null deviance: 21.60252
cat("Residual deviance:", model$deviance, "\n")
## Residual deviance: 9.32724
cat("AIC:", AIC(model), "\n")
## AIC: 41.35617

Residual deviance lebih kecil dari null deviance → model cukup baik.

  • Visualisasi Model
library(ggplot2)

ggplot(data, aes(x = jam, y = pelanggan)) +
  geom_point(color = "blue", size = 3) +
  geom_smooth(method = "glm", method.args = list(family = poisson(link = "log")), se = FALSE, color = "red") +
  labs(title = "Model GLM Poisson: Pelanggan vs Jam", x = "Jam", y = "Jumlah Pelanggan") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

  • Prediksi

Prediksi jumlah pelanggan pada jam ke-12:

jam_baru <- data.frame(jam = 12)
predict(model, newdata = jam_baru, type = "response")
##        1 
## 12.38555

Kesimpulan

  • GLM Poisson cocok untuk data jumlah kejadian (count).
  • Variabel waktu (jam) berpengaruh signifikan.
  • Setiap jam tambahan meningkatkan jumlah pelanggan sekitar 23%.
  • Model dapat digunakan untuk prediksi di masa depan.

Model Regresi Logistik

Regresi logistik adalah teknik statistik yang digunakan untuk memodelkan probabilitas dari suatu kejadian biner, yaitu hasil dengan dua kemungkinan, seperti ya/tidak atau 1/0. Berbeda dari regresi linear yang menghasilkan nilai kontinu, regresi logistik membatasi prediksi dalam rentang 0 hingga 1 menggunakan fungsi sigmoid.

Fungsi Sigmoid

Fungsi sigmoid digunakan untuk mengubah nilai input menjadi nilai probabilitas antara 0 dan 1. Rumus matematisnya adalah:

\[ \sigma(z) = \frac{1}{1 + e^{-z}} \]

Dalam regresi logistik, nilai \(z\) adalah kombinasi linear dari variabel prediktor:

\[ z = \beta_0 + \beta_1X_1 + \beta_2X_2 + \ldots + \beta_pX_p \]

Interpretasi Model

Model regresi logistik mengestimasi probabilitas keanggotaan dalam kelas tertentu. Jika hasil dari fungsi sigmoid lebih besar dari ambang batas (biasanya 0.5), maka model memprediksi kelas 1; jika lebih kecil, maka memprediksi kelas 0.

Contoh Kasus: Prediksi Kelayakan Kredit

Misalkan kita ingin memprediksi apakah seorang nasabah layak diberi kredit berdasarkan usia dan pendapatannya.

# Data simulasi
set.seed(123)
data_kredit <- data.frame(
  usia = sample(20:60, 100, replace = TRUE),
  pendapatan = sample(3:15, 100, replace = TRUE),
  layak_kredit = rbinom(100, 1, 0.5)
)

# Melatih model regresi logistik
model_logit <- glm(layak_kredit ~ usia + pendapatan, data = data_kredit, family = binomial())

# Ringkasan model
summary(model_logit)
## 
## Call:
## glm(formula = layak_kredit ~ usia + pendapatan, family = binomial(), 
##     data = data_kredit)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.72477    0.98411   0.736    0.461
## usia        -0.03074    0.01928  -1.594    0.111
## pendapatan   0.01672    0.06090   0.275    0.784
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 135.37  on 99  degrees of freedom
## Residual deviance: 132.73  on 97  degrees of freedom
## AIC: 138.73
## 
## Number of Fisher Scoring iterations: 4

Interpretasi Hasil

Deviasi Null (135.37): Deviasi ini menggambarkan seberapa baik model tanpa prediktor (hanya intercept) dalam menjelaskan data. Semakin rendah nilai deviasi ini, semakin baik model.

Deviasi Residual (132.73): Ini adalah deviasi setelah model dengan prediktor diterapkan. Model kita sedikit lebih baik daripada model null, meskipun penurunan deviasi relatif kecil.

AIC (138.73): AIC digunakan untuk membandingkan model yang berbeda. Model dengan AIC yang lebih rendah dianggap lebih baik, tetapi kita perlu membandingkannya dengan model lain untuk penilaian lebih lanjut.

Kesimpulan

Berdasarkan p-value yang lebih besar dari 0.05 untuk kedua variabel (usia dan pendapatan), kita gagal menolak hipotesis nol bahwa koefisien variabel-variabel ini adalah nol (tidak signifikan). Ini berarti bahwa dalam data simulasi ini, usia dan pendapatan tidak memiliki pengaruh yang signifikan terhadap kelayakan kredit pada tingkat signifikansi 0.05.

  • Visualisasi Fungsi Sigmoid

Berikut contoh kurva sigmoid:

sigmoid <- function(z) {
  1 / (1 + exp(-z))
}

curve(sigmoid(x), from = -10, to = 10, col = "blue", ylab = "Probabilitas", main = "Kurva Fungsi Sigmoid")

Keunggulan Regresi Logistik

  • Implementasi sederhana dan efisien untuk klasifikasi biner.
  • Cocok untuk data yang dapat dipisahkan secara linear.
  • Koefisien mudah diinterpretasi sebagai pengaruh log-odds.
  • Cocok digunakan dalam berbagai bidang: kesehatan, keuangan, pemasaran, dan lainnya.

Asumsi Model

  • Tidak ada multikolinearitas antar variabel prediktor.
  • Observasi bersifat independen.
  • Hubungan antara prediktor dan log-odds bersifat linear.

Model regresi Poisson

Regresi Poisson digunakan ketika variabel respons berupa data cacah (count data), yaitu bilangan bulat non-negatif seperti jumlah pelanggan, jumlah kecelakaan, jumlah kejadian, dan sebagainya.

Model ini merupakan bagian dari Generalized Linear Model (GLM) dengan asumsi bahwa distribusi dari variabel respons mengikuti distribusi Poisson.

Distribusi Poisson

Distribusi Poisson memiliki fungsi probabilitas:

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

Kita dapat menuliskan bentuk ini dalam format exponential family:

\[ 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) \]

Model regresi menjadi:

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

Fungsi inverse-nya adalah:

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

Estimasi Parameter

Estimasi parameter \(\beta\) dilakukan menggunakan metode Maximum Likelihood Estimation (MLE).

Fungsi log-likelihood untuk model Poisson:

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

Estimasi \(\beta\) diperoleh melalui metode numerik, misalnya iterasi Newton-Raphson.

**Contoh Kasus:* Jumlah Transaksi per Hari

Misalkan kita memiliki data simulasi mengenai jumlah transaksi yang terjadi per hari, tergantung pada nilai promosi harian yang diterapkan (nilai diskon). Kita ingin melihat bagaimana nilai promosi mempengaruhi jumlah transaksi.

set.seed(42)
n <- 200
promosi <- rnorm(n)
lambda <- exp(0.3 + 0.6 * promosi)
transaksi <- rpois(n, lambda)
data <- data.frame(transaksi, promosi)
head(data)
  • Estimasi Regresi Poisson
poisson_model <- glm(transaksi ~ promosi, data = data, family = poisson)
summary(poisson_model)
## 
## Call:
## glm(formula = transaksi ~ promosi, family = poisson, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.21817    0.06712   3.250  0.00115 ** 
## promosi      0.58748    0.06288   9.343  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 333.72  on 199  degrees of freedom
## Residual deviance: 244.46  on 198  degrees of freedom
## AIC: 587.36
## 
## Number of Fisher Scoring iterations: 5
  • Visualisasi Model
plot(data$promosi, data$transaksi, pch = 16, col = "darkgray", 
     main = "Data dan Hasil Prediksi Jumlah Transaksi")
newdata <- data.frame(promosi = seq(min(promosi), max(promosi), length.out = 100))
pred <- predict(poisson_model, newdata = newdata, type = "response")
lines(newdata$promosi, pred, col = "blue", lwd = 2)

  • Interpretasi Koefisien
exp(coef(poisson_model))
## (Intercept)     promosi 
##    1.243799    1.799442

Koefisien eksponensial dari variabel promosi menunjukkan faktor pengali terhadap rata-rata jumlah transaksi untuk setiap kenaikan satu unit promosi.

  • Visualisasi dengan ggplot2 (Kasus warpbreaks)
library(ggplot2)
poisson_model2 <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)
warpbreaks$predicted <- predict(poisson_model2, type = "response")

ggplot(warpbreaks, aes(x = tension, y = breaks, color = wool)) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  geom_point(aes(y = predicted), shape = 18, size = 3, color = "black") +
  facet_wrap(~wool) +
  labs(title = "Prediksi Jumlah Patahan Benang berdasarkan Wol dan Ketegangan",
       x = "Tingkat Ketegangan",
       y = "Jumlah Patahan Benang",
       color = "Jenis Wol") +
  theme_minimal()

Kesimpulan

  • Regresi Poisson sesuai untuk memodelkan data cacah.
  • Contoh simulasi menunjukkan hubungan eksponensial antara tingkat promosi dan jumlah transaksi.
  • Hasil model dapat diinterpretasikan dengan fungsi link logaritma dan eksponensial terhadap koefisien.
  • Visualisasi memberikan gambaran model yang baik dalam menggambarkan tren data.

Inferensi GLM

Dalam Generalized Linear Model (GLM), inferensi statistik berperan penting dalam menilai ketepatan dan reliabilitas model. Inferensi ini mencakup ekspektasi dan varians dari estimator parameter, serta distribusi asimptotik dari estimator tersebut.

Ekspektasi Estimator

Ekspektasi menunjukkan apakah suatu estimator tak bias, yaitu:

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

Dalam GLM, estimator maksimum likelihood (MLE) dari \(\hat{\beta}\) bersifat asymptotically unbiased, artinya ketika ukuran sampel cukup besar, rata-rata nilai \(\hat{\beta}\) mendekati nilai parameter sebenarnya \(\beta\).

Varians Estimator

Varians menunjukkan presisi dari estimator parameter. Varians dari \(\hat{\beta}\) dalam GLM mendekati:

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

di mana \(\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 normal asimptotik ini menjadi dasar bagi:

Varians dalam GLM Tidak Konstan

Berbeda dengan regresi linear (OLS) yang mengasumsikan:

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

Dalam GLM:

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

Contoh:

Contoh Simulasi dan Estimasi Varians Model

set.seed(456)
x <- rnorm(100)
p <- exp(0.5 + 0.8 * x) / (1 + exp(0.5 + 0.8 * x))  # Probabilitas antara 0 dan 1
y <- rbinom(100, size = 1, prob = p)  # Data binomial (sukses/gagal)
data <- data.frame(x, y)

# Model GLM dengan distribusi Binomial
model <- glm(y ~ x, family = binomial, data = data)
summary(model)
## 
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.5157     0.2234   2.308 0.020984 *  
## x             0.8652     0.2525   3.427 0.000611 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 131.79  on 99  degrees of freedom
## Residual deviance: 117.56  on 98  degrees of freedom
## AIC: 121.56
## 
## Number of Fisher Scoring iterations: 3

Model ini memperkirakan hubungan antara prediktor \(x\) dan data cacah \(y\) menggunakan regresi Poisson. Koefisien model dapat diuji dengan uji Wald menggunakan estimasi standar error.

# Ekstrak koefisien dan varians kovarian
coef(model)
## (Intercept)           x 
##   0.5157375   0.8652363
vcov(model)
##             (Intercept)           x
## (Intercept) 0.049920864 0.004197341
## x           0.004197341 0.063750699

Kovarian dan standar error ini digunakan untuk membangun confidence interval dan menghitung statistik Wald.

Visualisasi

Mencari Ekspektasi dan Varians dalam GLM

Ekspektasi

Dalam konteks Generalized Linear Model (GLM), ekspektasi dari variabel respon \(Y\) dapat diperoleh dari fungsi momen. Ekspektasi ini menyatakan nilai tengah atau rata-rata dari distribusi \(Y\).

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

Distribusi dalam GLM umumnya berasal dari keluarga eksponensial, yang bentuk umumnya adalah:

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

Namun, dalam bentuk yang lebih sederhana sering dituliskan sebagai:

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

Turunan pertama dari log-likelihood terhadap parameter \(\theta\) disebut fungsi skor:

\[ U(\theta) = \frac{\partial \ell}{\partial \theta} = y - b'(\theta) \]

Ekspektasi dari Turunan Pertama

Karena \(E[U(\theta)] = 0\), maka:

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

Sehingga diperoleh bahwa:

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

Varians

Untuk mencari varians dari \(Y\), kita turunkan kembali fungsi log-likelihood terhadap \(\theta\), menghasilkan turunan kedua:

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

Karena \(Var(Y) = E[(Y - \mu)^2]\), maka berdasarkan properti distribusi eksponensial:

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

Di mana: - \(\phi\) adalah dispersi, - \(V(\mu)\) adalah fungsi varians yang tergantung pada \(\mu\) (berbeda-beda tergantung distribusinya, misalnya \(V(\mu) = \mu\) untuk Poisson, atau \(V(\mu) = \mu(1 - \mu)\) untuk binomial).

Contoh Kasus

Misalkan kita memiliki data jumlah pelanggan yang datang ke sebuah kedai kopi setiap jam, dan jumlah tersebut mengikuti distribusi Poisson, yang termasuk dalam keluarga eksponensial.

Distribusi Poisson:

\[ f(y; \lambda) = \frac{e^{-\lambda} \lambda^y}{y!} \]

Jika kita nyatakan dalam bentuk eksponensial:

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

Kita bisa mengenali bahwa: - \(\theta = \log \lambda\) - \(b(\theta) = e^\theta\) - \(\mu = b'(\theta) = e^\theta = \lambda\)

Dan:

\[ Var(Y) = b''(\theta) = e^\theta = \lambda = \mu \]

# Simulasi data Poisson
set.seed(123)
y <- rpois(100, lambda = 4)

# Ekspektasi dan varians dari sampel
mean(y)  # estimasi ekspektasi
## [1] 4.09
var(y)   # estimasi varians
## [1] 3.820101

*Kesimpulan**

  • Distribusi Poisson adalah salah satu anggota dari keluarga eksponensial, sehingga cocok digunakan dalam GLM.

  • Simulasi dengan 100 observasi menghasilkan estimasi:

Ekspektasi mendekati 4 (tepatnya r round(mean_y, 2))

Varians mendekati 4 (tepatnya r round(var_y, 2))

  • Perbedaan kecil disebabkan oleh fluktuasi acak dari proses simulasi. Jika ukuran sampel diperbesar, nilai estimasi akan semakin mendekati nilai sebenarnya.

Diagnostik Model GLM

Diagnostik dalam GLM bertujuan untuk mengevaluasi apakah model yang telah dibangun sudah sesuai dengan data. Evaluasi ini dilakukan melalui:

  • Uji formal menggunakan statistik devians atau Chi-kuadrat.
  • Visualisasi seperti grafik antara nilai prediksi dan nilai aktual (observasi).

Statistik Devians

Statistik devians digunakan untuk mengukur seberapa baik model yang dibangun jika dibandingkan dengan model yang paling cocok (saturated model). Rumus umum devians adalah:

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

Interpretasi: - Nilai devian besar → model kurang cocok. - Nilai devian kecil → model lebih cocok dengan data. - Biasanya dibandingkan dengan model alternatif, atau digunakan sebagai dasar dalam uji devians.

Statistik Chi-Kuadrat Pearson

Statistik ini mengukur sejauh mana model yang dibangun lebih baik dibandingkan dengan tidak adanya model (model nol). Rumus statistik Chi-Square Pearson:

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

Interpretasi: - Nilai \(X^2\) yang besar dapat menunjukkan bahwa model kurang cocok. - Jika hasil uji signifikan → model memberikan penyesuaian yang lebih baik daripada model nol.

Catatan Penting

  • Untuk data berkelompok, devians dan statistik Pearson mengikuti distribusi Chi-kuadrat.
  • Untuk data tidak berkelompok (data individu), distribusinya tidak dapat dipastikan.
  • Devians diminimalkan oleh Maximum Likelihood Estimation (MLE), sehingga cocok digunakan sebagai alat evaluasi model.

Analisis Residual

  • Residual adalah selisih antara nilai yang diobservasi dengan nilai yang diprediksi oleh model.
  • Digunakan untuk mengevaluasi fit model dan untuk mengidentifikasi observasi yang mencurigakan (outliers).
  • Residual juga bisa diplot untuk melihat apakah terdapat pola sistematis.

Contoh Kasus: Model GLM Poisson

Kita simulasikan data jumlah pengunjung toko tiap jam (mengikuti distribusi Poisson), lalu evaluasi model menggunakan statistik devians, Chi-kuadrat, dan residual.

set.seed(123)

# Simulasi data
jam <- 1:100
pengunjung <- rpois(100, lambda = 5 + 0.1 * jam)  # model dengan tren meningkat
data <- data.frame(jam, pengunjung)

# Bangun model GLM Poisson
model <- glm(pengunjung ~ jam, family = poisson(link = "log"), data = data)

# Ringkasan model
summary(model)
## 
## Call:
## glm(formula = pengunjung ~ jam, family = poisson(link = "log"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.830677   0.071980  25.433  < 2e-16 ***
## jam         0.008476   0.001123   7.548 4.43e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 148.835  on 99  degrees of freedom
## Residual deviance:  90.846  on 98  degrees of freedom
## AIC: 501.2
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

  • Intercept: 1.83 → rata-rata log jumlah pengunjung saat jam = 0.

  • Koefisien jam: 0.0085 → menunjukkan bahwa jumlah pengunjung meningkat seiring waktu, karena nilai positif dan signifikan.

Kedua koefisien signifikan secara statistik dengan nilai-p yang sangat kecil, artinya ada hubungan yang cukup kuat antara prediktor jam dengan jumlah pengunjung.

  • Evaluasi Model
# Statistik devians
deviance(model)
## [1] 90.84634
# Statistik Pearson
pearson_resid <- residuals(model, type = "pearson")
chi_sq_pearson <- sum(pearson_resid^2)
chi_sq_pearson
## [1] 89.44523

Interpretasi

  • Devians

Dibandingkan dengan null deviance (148.83), model berhasil menurunkan devians cukup signifikan, menunjukkan model lebih baik dari tanpa prediktor.

Nilai ini dapat diuji secara formal dengan uji Likelihood Ratio, tapi secara umum, penurunan devians ini menunjukkan fit yang lebih baik.

  • Chi-Kuadrat Pearson

Nilai ini mirip dengan devians, yang menunjukkan tidak ada masalah besar dalam fit model.

Untuk Poisson dengan n = 100, nilai devians dan Pearson Chi-Square mendekati derajat kebebasan (98) → ini pertanda model sudah cocok dengan baik.

  • Visualisasi Residual
# Plot residual Pearson vs fitted values
plot(fitted(model), pearson_resid,
     xlab = "Fitted values", ylab = "Pearson Residuals",
     main = "Plot Pearson Residuals")
abline(h = 0, col = "red")

Kesimpulan

  • Model Poisson yang dibangun dengan prediktor jam sudah cukup baik dan cocok berdasarkan statistik devians dan Pearson.

  • Koefisien signifikan, artinya jam berkontribusi terhadap jumlah pengunjung.

  • Plot residual menyebar acak, yang mengindikasikan tidak ada pola yang terlewatkan model.

Secara keseluruhan, tidak ada bukti kuat bahwa model tidak cocok, sehingga dapat digunakan untuk interpretasi atau prediksi lebih lanjut.

Detail Metode Estimasi dan Inferensi Regresi Logistik

Regresi logistik digunakan untuk memodelkan probabilitas dari suatu peristiwa (biasanya dikodekan sebagai 1) berdasarkan satu atau lebih variabel prediktor. Model ini sangat berguna ketika variabel respon bersifat biner, seperti “ya/tidak”, “sukses/gagal”, atau “positif/negatif”. Estimasi parameter dilakukan menggunakan metode Maximum Likelihood Estimation (MLE), karena model logistik bersifat non-linear dalam parameter.

Fungsi Model Logistik

Fungsi model logistik menyatakan probabilitas \(\pi(x)\) dari keberhasilan (misalnya \(y=1\)) sebagai:

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

Fungsi 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

Newton-Raphson digunakan untuk mencari nilai parameter \(\beta\) yang memaksimumkan log-likelihood.

Turunan Pertama (Score Function):

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

Turunan Kedua (Hessian Matrix):

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

Iterasi Newton-Raphson:

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

Contoh Kasus: Prediksi Kelayakan Kredit Misalnya kita ingin memprediksi apakah seseorang akan membayar pinjamannya berdasarkan skor kredit. Data disimulasikan sebagai berikut:

set.seed(1)
score <- rnorm(100, mean = 600, sd = 50)
beta_true <- c(-8, 0.01)
X <- cbind(1, score)
eta <- X %*% beta_true
p <- 1 / (1 + exp(-eta))
y <- rbinom(100, 1, p)

# Estimasi Newton-Raphson manual
beta <- c(0, 0)
tol <- 1e-6
max_iter <- 100
for (i in 1:max_iter) {
  eta <- X %*% beta
  pi_hat <- 1 / (1 + exp(-eta))
  W <- diag(as.numeric(pi_hat * (1 - pi_hat)))
  z <- eta + solve(W) %*% (y - pi_hat)
  beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
  if (sum(abs(beta_new - beta)) < tol) break
  beta <- beta_new
}
beta
##               [,1]
##       -12.08946093
## score   0.01636156

Interpretasi

  • Intercept yang negatif besar menunjukkan bahwa individu dengan skor kredit sangat rendah memiliki probabilitas keberhasilan (membayar pinjaman) yang sangat kecil.
  • Koefisien skor kredit sebesar 0.0164 menunjukkan bahwa setiap peningkatan 1 poin skor kredit akan meningkatkan log-odds keberhasilan sebesar 0.0164. Dalam konteks probabilitas, ini berarti semakin tinggi skor kredit, semakin besar peluang seseorang untuk membayar pinjamannya.
  • Nilai koefisien yang diperoleh relatif dekat dengan nilai parameter sebenarnya yang digunakan saat simulasi data (β₀ = -8 dan β₁ = 0.01), yang menunjukkan bahwa algoritma estimasi Newton-Raphson bekerja dengan baik.

Inferensi Parameter

Uji Wald

Digunakan untuk mengetahui apakah suatu parameter \(\beta_j\) secara statistik berbeda dari nol.

Hipotesis: - \(H_0\): \(\beta_j = 0\) - \(H_1\): \(\beta_j \ne 0\)

Jika estimator \(\hat{\beta}_j\) berdistribusi normal, maka statistik uji adalah:

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

Statistik Wald dihitung sebagai:

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

Contoh Kasus: Prediksi Respon Terhadap Promo

Misalkan kita ingin tahu apakah intensitas iklan (jumlah jam seseorang melihat iklan) memengaruhi keputusan membeli produk (1 = beli, 0 = tidak beli).

set.seed(123)
n <- 100
iklan <- runif(n, 0, 10)
log_odds <- -1 + 0.6 * iklan
p <- 1 / (1 + exp(-log_odds))
beli <- rbinom(n, 1, p)
data <- data.frame(iklan, beli)

model <- glm(beli ~ iklan, data = data, family = binomial)
summary(model)
## 
## Call:
## glm(formula = beli ~ iklan, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.4405     0.4940  -0.892 0.372546    
## iklan         0.4834     0.1326   3.645 0.000268 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 97.245  on 99  degrees of freedom
## Residual deviance: 77.838  on 98  degrees of freedom
## AIC: 81.838
## 
## Number of Fisher Scoring iterations: 5
# Langkah-langkah Uji Wald
beta_hat <- coef(model)["iklan"]
se_beta <- summary(model)$coefficients["iklan", "Std. Error"]
Z <- beta_hat / se_beta
Wald_stat <- Z^2
p_value <- 1 - pchisq(Wald_stat, df = 1)

Z
##    iklan 
## 3.644741
Wald_stat
##    iklan 
## 13.28414
p_value
##        iklan 
## 0.0002676614

Interpretasi:

  • Koefisien 0.4834 untuk variabel iklan menunjukkan bahwa setiap tambahan satu jam melihat iklan meningkatkan log-odds seseorang membeli produk sebesar 0.4834.
  • Dalam bentuk odds, nilai ini ekuivalen dengan peningkatan odds sebesar \(\exp(0.4834) \approx 1.621\), artinya setiap tambahan 1 jam melihat iklan meningkatkan peluang pembelian sebesar 62.1%.
  • Nilai p untuk iklan adalah 0.000268, yang jauh lebih kecil dari 0.05, sehingga kita tolak hipotesis nol dan menyimpulkan bahwa intensitas iklan berpengaruh signifikan terhadap keputusan pembelian.
  • Statistik Wald sebesar 13.28 dengan p-value 0.00027 menguatkan kesimpulan tersebut.

Uji Likelihood Ratio (Chi-Square)

Digunakan untuk membandingkan model penuh (dengan prediktor) dan model null (tanpa prediktor).

# Model null
model_null <- glm(beli ~ 1, data = data, family = binomial)

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

Interpretasi:

  • Deviance model null = 97.245, deviance model penuh = 77.838, selisih deviance sebesar 19.406.
  • Dengan derajat bebas 1, p-value yang dihasilkan < 0.001, sehingga kita tolak hipotesis nol bahwa model null cukup baik.
  • Artinya, penambahan variabel iklan secara signifikan meningkatkan kualitas model.

Evaluasi Kebaikan Model

  1. Akaike Information Criterion (AIC) AIC menghitung trade-off antara goodness-of-fit dan kompleksitas model. Semakin kecil AIC, semakin baik model.
AIC(model)
## [1] 81.83842
  1. Bayesian Information Criterion (BIC) BIC mirip dengan AIC, tetapi memberikan penalti lebih besar terhadap jumlah parameter:
BIC(model)
## [1] 87.04876

Interpretasi:

  • Nilai AIC model = 81.84 dan BIC = 87.05.
  • Nilai-nilai ini bisa digunakan untuk membandingkan dengan model alternatif. Semakin kecil nilainya, semakin baik model tersebut dari segi keseimbangan antara fit dan kompleksitas.

Kesimpulan:

  • Estimasi parameter regresi logistik dilakukan melalui pendekatan MLE menggunakan metode iteratif Newton-Raphson.
  • Inferensi parameter dapat dilakukan dengan uji Wald dan uji Likelihood Ratio (Chi-Square).
  • Evaluasi performa model bisa dilakukan dengan AIC dan BIC untuk memilih model terbaik dengan kompleksitas minimum.
  • Contoh-contoh yang digunakan menunjukkan bagaimana regresi logistik dapat diaplikasikan pada kasus riil seperti prediksi kelayakan kredit dan efektivitas iklan, lengkap dengan interpretasi hasil estimasi dan pengujian statistik.

Detail Metode Estimasi dan Inferensi Regresi Logistik

Model regresi Poisson digunakan untuk memodelkan data count (jumlah kejadian) yang mengikuti distribusi Poisson. Model ini sangat berguna untuk data diskrit non-negatif, seperti jumlah pelanggan, jumlah kecelakaan, dan lainnya.

Estimasi parameter dilakukan menggunakan Maximum Likelihood Estimation (MLE), dan inferensi dilakukan menggunakan uji Wald serta Likelihood Ratio Test. Salah satu pendekatan numerik yang digunakan untuk mendapatkan estimasi parameter adalah Iteratively Reweighted Least Squares (IRLS).

Model Regresi Poisson

Distribusi Poisson didefinisikan sebagai:

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

Dengan \(\lambda_i = \mathbb{E}[Y_i]\). Model regresi Poisson menggunakan fungsi link logaritma:

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

Fungsi Log-Likelihood

Fungsi log-likelihood yang digunakan untuk estimasi parameter adalah:

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

Dengan \(\lambda_i = \exp(x_i^\top \beta)\).

Estimasi dengan Iteratively Reweighted Least Squares (IRLS)

Langkah-langkah IRLS sebagai berikut:

Langkah 1: Definisikan model regresi

\[ \eta_i = x_i^\top \beta \Rightarrow \lambda_i = \exp(\eta_i) \]

Langkah 2: Bentuk log-likelihood

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

Langkah 3: Rumus iterasi IRLS

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

Dengan:

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

Contoh Kasus: Estimasi Manual (IRLS)

# Simulasi data
set.seed(456)
n <- 100
x <- rnorm(n, mean = 2, sd = 1)
X <- cbind(1, x)  # Tambah intercept
beta_true <- c(0.7, 1.2)
eta <- X %*% beta_true
lambda <- exp(eta)
y <- rpois(n, lambda)
  • IRLS Manual Step-by-Step
beta <- c(0, 0)
tol <- 1e-6
max_iter <- 100

for (i in 1:max_iter) {
  eta <- X %*% beta
  lambda <- exp(eta)
  W <- diag(as.numeric(lambda))
  z <- eta + (y - lambda) / lambda
  beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z

  if (sum(abs(beta_new - beta)) < tol) {
    cat("Konvergen pada iterasi ke-", i, "\n")
    break
  }

  beta <- beta_new
}

beta  # Estimasi parameter akhir
##         [,1]
##   -163.22640
## x   53.65828

Interpretasi

Didapatkan hasil nilai intercept sangat negatif dan slope sangat tinggi, menandakan kemungkinan numerik yang tidak stabil atau outlier dalam data. Oleh karena itu, estimasi dari glm() lebih dipercaya untuk interpretasi.

  • Estimasi dengan glm()
model_glm <- glm(y ~ x, family = poisson)
summary(model_glm)
## 
## Call:
## glm(formula = y ~ x, family = poisson)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.68664    0.05658   12.14   <2e-16 ***
## x            1.20663    0.01726   69.90   <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: 6109.69  on 99  degrees of freedom
## Residual deviance:  103.49  on 98  degrees of freedom
## AIC: 610.73
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

Model regresi Poisson yang dipasang menggunakan fungsi glm() menghasilkan:

  • Intercept: 0.68664
  • Koefisien x: 1.20663

Artinya, setiap kenaikan satu unit pada variabel prediktor x, log-rata-rata kejadian (jumlah y) meningkat sebesar 1.20663. Jika diubah ke skala asli (eksponensial), maka jumlah kejadian meningkat sekitar:

\[ \exp(1.20663) \approx 3.34 \text{ kali} \]

  • Uji Hipotesis
  1. Uji Wald

Untuk menguji hipotesis nol \(H_0: \beta_1 = 0\):

coef_val <- coef(model_glm)[2]
se_val <- summary(model_glm)$coefficients[2, 2]
wald_z <- coef_val / se_val
wald_chisq <- wald_z^2
p_value <- 1 - pchisq(wald_chisq, df = 1)

cat("Z:", wald_z, "\nChi-Square:", wald_chisq, "\np-value:", p_value)
## Z: 69.8977 
## Chi-Square: 4885.688 
## p-value: 0

Interpretasi

Uji Wald terhadap koefisien x menghasilkan statistik Z sebesar 69.90 dengan p-value ≈ 0. Hasil ini menunjukkan bahwa koefisien x signifikan secara statistik, sehingga menolak hipotesis nol \(H_0: \beta_1 = 0\).

  1. Uji Likelihood Ratio
model_null <- glm(y ~ 1, family = poisson)
anova(model_null, model_glm, test = "Chisq")

Interpretasi

Uji Likelihood Ratio membandingkan model penuh (y ~ x) dengan model null (y ~ 1):

  • Nilai deviance berkurang dari 6109.7 menjadi 103.5
  • Perbedaan deviance sebesar 6006.2 sangat signifikan (p < 2.2e-16)

Hal ini menunjukkan bahwa variabel x secara signifikan meningkatkan kemampuan model dalam menjelaskan variasi data.

  • Evaluasi Model

AIC dan BIC

AIC(model_glm)
## [1] 610.7341
BIC(model_glm)
## [1] 615.9445

Interpretasi

Nilai AIC dan BIC relatif rendah, menunjukkan model cukup baik dalam menjelaskan data dengan kompleksitas yang wajar.

Kesimpulan

Model regresi Poisson dapat diestimasi baik secara manual (IRLS) maupun menggunakan fungsi glm(). Dalam praktiknya, glm() lebih stabil dan akurat. Berdasarkan hasil simulasi:

  • Koefisien prediktor x signifikan secara statistik
  • Model dengan x jauh lebih baik dari model tanpa prediktor (berdasarkan deviance dan uji LR)
  • Model memiliki AIC dan BIC yang kompetitif

Regresi Poisson sangat cocok untuk memodelkan data count seperti ini.

Contoh Kasus : Analisis Conditional Independence antar Produk pada Data Market Basket Transaksi Groceries

Pendahuluan

Analisis keranjang belanja (market basket analysis) adalah pendekatan eksploratif dalam data mining untuk menemukan pola hubungan antar produk dalam transaksi pelanggan. Salah satu pendekatan lanjutannya adalah mengevaluasi kemandirian bersyarat atau conditional independence, yaitu apakah dua produk tetap berhubungan setelah mempertimbangkan pengaruh produk ketiga.

Urgensi

Memahami keterkaitan bersyarat antar produk bisa membantu:

  • Penataan produk yang efisien
  • Perencanaan promosi kombinasi produk
  • Segmentasi pelanggan berdasarkan pola konsumsi

Dataset

Dataset yang digunakan adalah “Groceries”, berisi 9835 transaksi belanja dengan ratusan produk unik. Setiap baris adalah satu transaksi, dan setiap transaksi bisa berisi satu atau lebih produk. Dataset bersumber dari website Kaggle.com

library(arules)
## Loading required package: Matrix
## 
## Attaching package: 'arules'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following objects are masked from 'package:base':
## 
##     abbreviate, write
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.1
## ✔ readr     2.1.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()     masks Matrix::expand()
## ✖ tidyr::extract()    masks magrittr::extract()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::group_rows() masks kableExtra::group_rows()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ tidyr::pack()       masks Matrix::pack()
## ✖ arules::recode()    masks dplyr::recode()
## ✖ purrr::set_names()  masks magrittr::set_names()
## ✖ tidyr::unpack()     masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

load data

# Load libraries
library(arules)
library(tidyverse)

# Baca data transaksi
data <- read.transactions(
  file = "C:/Users/asus/Downloads/groceries - groceries.csv",
  format = "basket",
  sep = ","
)
## Warning in readLines(file, encoding = encoding): incomplete final line found on
## 'C:/Users/asus/Downloads/groceries - groceries.csv'
# Ringkasan data
summary(data)
## transactions as itemMatrix in sparse format with
##  9836 rows (elements/itemsets/transactions) and
##  231 columns (items) and a density of 0.0234297 
## 
## most frequent items:
##       whole milk                1 other vegetables       rolls/buns 
##             2513             2159             1903             1809 
##             soda          (Other) 
##             1715            43136 
## 
## element (itemset/transaction) length distribution:
## sizes
##    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17 
## 2159 1643 1299 1005  855  645  545  438  350  246  182  117   78   77   55   46 
##   18   19   20   21   22   23   24   25   27   28   29   30   33 
##   29   14   14    9   11    4    6    1    1    1    1    3    2 
## 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.000   3.000   4.000   5.412   7.000  33.000 
## 
## includes extended item information - examples:
##   labels
## 1      1
## 2     10
## 3     11
# Visualisasi item terbanyak
itemFrequencyPlot(data, topN = 10, type = "absolute", main = "10 Produk Paling Sering Dibeli")

Interpretasi Grafik Produk Populer

Dari grafik batang di atas, terlihat bahwa produk yang paling sering dibeli adalah: - whole milk, - other vegetables, - rolls/buns, - soda, - dan yogurt, di antara item-item lainnya.

Produk-produk ini menarik untuk dianalisis karena frekuensinya tinggi, yang berarti kita punya cukup data untuk mengevaluasi hubungan di antaranya.

Pemilihan Variabel untuk Analisis Conditional Independence

Untuk melakukan analisis kemandirian bersyarat, kita pilih tiga item berdasarkan popularitas dan potensi keterkaitan logis antar produk:

  • X: rolls/buns (karbohidrat/pastry)
  • Y: soda (minuman ringan)
  • Z: other vegetables (sayuran lain, sebagai indikator pola makan sehat)

Alasan Pemilihan:

  • Kita ingin tahu apakah pelanggan yang membeli roti (rolls/buns) juga cenderung membeli soda, dan apakah hubungan ini masih berlaku setelah dikendalikan oleh keberadaan sayur (other vegetables).
  • Ini mencerminkan analisis pola konsumsi sehat vs tidak sehat, yang relevan untuk pemasaran atau segmentasi konsumen.

Dengan pendekatan ini, kita bisa menguji apakah preferensi produk terjadi karena kebiasaan konsisten atau dipengaruhi oleh faktor lain (Z).

  1. Analisis Conditional Independence
# Ubah ke data frame biner
binary_df <- as(data, "matrix") %>% as.data.frame()

# Definisikan item
X <- binary_df$`rolls/buns`
Y <- binary_df$soda
Z <- binary_df$`other vegetables`

# Bagi data berdasarkan Z
Z_true <- binary_df[Z == 1, ]
Z_false <- binary_df[Z == 0, ]

# Buat tabel kontingensi saat Z = 1
tab_Z_true <- table(Z_true$`rolls/buns`, Z_true$soda)

# Buat tabel kontingensi saat Z = 0
tab_Z_false <- table(Z_false$`rolls/buns`, Z_false$soda)

# Uji chi-square
chi_true <- chisq.test(tab_Z_true)
chi_false <- chisq.test(tab_Z_false)

chi_true
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab_Z_true
## X-squared = 14.271, df = 1, p-value = 0.0001583
chi_false
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab_Z_false
## X-squared = 7.5603, df = 1, p-value = 0.005967

Interpretasi

Hasil uji Chi-Square memberikan dua skenario berdasarkan variabel kontrol other vegetables :

  • Saat Z = 1 (konsumen membeli sayur):
    • Nilai X-squared = 14.271 dengan p-value = 0.0001583
    • Ini menunjukkan ada hubungan yang signifikan antara pembelian rolls/buns dan soda ketika konsumen juga membeli other vegetables.
  • Saat Z = 0 (konsumen tidak membeli sayur):
    • Nilai X-squared = 7.5466 dengan p-value = 0.006012
    • Hubungan antara rolls/buns dan soda juga signifikan, meskipun kekuatannya lebih rendah dibandingkan saat Z = 1.

Kesimpulan Karena kedua kondisi (Z = 1 dan Z = 0) sama-sama menghasilkan p-value < 0.05, maka dapat disimpulkan bahwa:

Pembelian rolls/buns dan soda TIDAK independen secara kondisional terhadap other vegetables.

Artinya, hubungan antara dua produk ini tetap signifikan bahkan setelah mempertimbangkan apakah pelanggan membeli sayuran atau tidak. Dengan kata lain, ada asosiasi kuat yang tidak tereliminasi oleh faktor kontrol tersebut.

Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio 1

Simulasi Data

Data terdiri dari 500 observasi yang merepresentasikan individu dengan berbagai karakteristik perjalanan. Variabel respons berupa keputusan membeli asuransi (ya atau tidak), sementara variabel prediktor terdiri atas campuran skala: nominal (jenis perjalanan: domestik atau internasional), ordinal (tingkat kepanikan terhadap risiko perjalanan: Low, Medium, High, Very High), dan rasio (durasi perjalanan dalam hari). Simulasi ini memungkinkan analisis hubungan antara karakteristik perjalanan dan keputusan pembelian asuransi dengan pendekatan statistik yang fleksibel dan dapat dikendalikan.

library(tibble)
set.seed(2025)
n <- 500

# Nominal: Jenis perjalanan
travel_type <- sample(c("Domestic", "International"), n, replace = TRUE)

# Ordinal: Tingkat kepanikan terhadap risiko
risk_level <- sample(c("Low", "Medium", "High", "Very High"), n, 
                     replace = TRUE, prob = c(0.3, 0.4, 0.2, 0.1))

# Rasio: Durasi perjalanan dalam hari
trip_duration <- rpois(n, lambda = 7) + 1  # durasi minimal 1 hari

# Konversi risk_level jadi nilai ordinal numerik
risk_numeric <- as.numeric(factor(risk_level, 
                         levels = c("Low", "Medium", "High", "Very High"), 
                         ordered = TRUE))

# Model logit
logit_p <- -2 + 0.6 * (travel_type == "International") + 
               0.5 * risk_numeric + 
               0.1 * trip_duration

# Probabilitas membeli asuransi
p <- 1 / (1 + exp(-logit_p))

# Respons: membeli atau tidak
buy_insurance <- rbinom(n, 1, p)

# Gabungkan ke dalam data frame
sim_data <- tibble::tibble(buy_insurance, travel_type, risk_level, trip_duration)
head(sim_data)

Eksplorasi Data

library(dplyr)

sim_data %>% 
  dplyr::group_by(buy_insurance) %>% 
  dplyr::summarise(Jumlah = n(), Rata2_Durasi = mean(trip_duration))

Interpretasi Jumlah individu yang membeli dan tidak membeli asuransi adalah sama (masing-masing 250). Rata-rata durasi perjalanan bagi yang membeli asuransi (8.436 hari) sedikit lebih panjang dibandingkan yang tidak membeli (7.584 hari), yang mengindikasikan bahwa semakin lama perjalanan, kemungkinan membeli asuransi cenderung meningkat.

Perlakuan Variabel Ordinal

Perlakuan Risk Level sebagai Nominal

sim_data_nominal <- sim_data %>%
  mutate(
    risk_level = factor(risk_level, levels = c("Low", "Medium", "High", "Very High"))
  )

model_nominal <- glm(buy_insurance ~ travel_type + risk_level + trip_duration, 
                     data = sim_data_nominal, family = binomial)
summary(model_nominal)
## 
## Call:
## glm(formula = buy_insurance ~ travel_type + risk_level + trip_duration, 
##     family = binomial, data = sim_data_nominal)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -1.77335    0.35021  -5.064 4.11e-07 ***
## travel_typeInternational  0.77484    0.19125   4.051 5.09e-05 ***
## risk_levelMedium          0.35464    0.22307   1.590  0.11187    
## risk_levelHigh            0.81989    0.27575   2.973  0.00295 ** 
## risk_levelVery High       1.81932    0.37303   4.877 1.08e-06 ***
## trip_duration             0.11289    0.03603   3.133  0.00173 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 693.15  on 499  degrees of freedom
## Residual deviance: 633.92  on 494  degrees of freedom
## AIC: 645.92
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

Intercept: Log-odds membeli asuransi untuk referensi (perjalanan domestik, risiko Low, durasi = 0), yaitu -1.773.

travel_typeInternational: Peluang membeli asuransi lebih tinggi secara signifikan untuk perjalanan internasional (OR ≈ exp(0.775) ≈ 2.17).

risk_level:

Medium: Tidak signifikan.

High: Signifikan (p < 0.01), peluang meningkat.

Very High: Sangat signifikan (p < 0.001), peluang membeli meningkat tajam.

trip_duration: Semakin lama perjalanan, peluang membeli asuransi meningkat signifikan.

Perlakuan Risk Level sebagai Numeric

sim_data_numeric <- sim_data %>%
  mutate(
    risk_numeric = case_when(
      risk_level == "Low" ~ 1,
      risk_level == "Medium" ~ 2,
      risk_level == "High" ~ 3,
      risk_level == "Very High" ~ 4
    )
  )

model_numeric <- glm(buy_insurance ~ travel_type + risk_numeric + trip_duration,
                     data = sim_data_numeric, family = binomial)
summary(model_numeric)
## 
## Call:
## glm(formula = buy_insurance ~ travel_type + risk_numeric + trip_duration, 
##     family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -2.3974     0.3873  -6.190 6.00e-10 ***
## travel_typeInternational   0.7671     0.1903   4.031 5.55e-05 ***
## risk_numeric               0.5292     0.1014   5.219 1.80e-07 ***
## trip_duration              0.1136     0.0359   3.165  0.00155 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 693.15  on 499  degrees of freedom
## Residual deviance: 635.97  on 496  degrees of freedom
## AIC: 643.97
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

Perlakuan ordinal menghasilkan koefisien yang signifikan (p < 0.001) → semakin tinggi tingkat risiko (Low sampai Very High), log-odds membeli asuransi meningkat secara linier.

Hal ini lebih cocok jika tingkatan risiko memang teratur secara logika dan pengaruhnya juga bertingkat.

Perbandingan Model

nullmod <- glm(buy_insurance ~ 1, data = sim_data, family = binomial)

# AIC
AIC(model_nominal)
## [1] 645.9207
AIC(model_numeric)
## [1] 643.9686

Interpretasi

Model dengan risk level sebagai numeric memiliki nilai AIC lebih rendah, artinya model ini lebih baik secara keseluruhan (lebih fit dan lebih parsimonious) dibandingkan pendekatan nominal.

Goodness-of-Fit

1 - (logLik(model_nominal)/logLik(nullmod))
## 'log Lik.' 0.0854457 (df=6)
1 - (logLik(model_numeric)/logLik(nullmod))
## 'log Lik.' 0.08249125 (df=4)

Interpretasi

Nilai ini menyerupai pseudo R². Artinya, baik model nominal maupun numeric hanya menjelaskan sekitar 8–9% variasi dalam keputusan membeli asuransi.

Ini umum pada model logistik karena variabilitas seringkali besar dan kompleks.

Visualisasi Prediksi

library(ggplot2)
library(dplyr)

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

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

Model nominal

ggplot(sim_data_nominal, aes(x = trip_duration, y = predicted, color = risk_level)) +
  geom_point(alpha = 0.5) +
  labs(title = "Prediksi Probabilitas Pembelian Asuransi (Risk Level sebagai Dummy)",
       x = "Durasi Perjalanan (hari)", y = "Probabilitas Membeli") +
  theme_minimal()

Interpretasi

Grafik menunjukkan bahwa semakin lama durasi perjalanan, probabilitas membeli asuransi meningkat untuk semua level risiko.

Perbedaan antara level risiko terlihat jelas: “Very High” menghasilkan probabilitas tertinggi di semua durasi.

Model numeric

ggplot(sim_data_numeric, aes(x = trip_duration, y = predicted, color = as.factor(risk_numeric))) +
  geom_point(alpha = 0.5) +
  labs(title = "Prediksi Probabilitas Pembelian Asuransi (Risk Level sebagai Skor)",
       x = "Durasi Perjalanan (hari)", y = "Probabilitas Membeli") +
  theme_minimal()

Interpretasi

Tren yang serupa seperti pada model nominal, tapi di sini level risiko diperlakukan sebagai skala bertingkat.

Prediksi cenderung lebih “halus” dan linier karena risk level dipetakan menjadi 1–4.

Ringkasan model nominal

library(broom)
library(knitr)
library(kableExtra)

tidy(model_nominal) %>%
  mutate(estimate = round(estimate, 3),
         std.error = round(std.error, 3),
         p.value = round(p.value, 3)) %>%
  kable(format = "html", caption = "Ringkasan Koefisien Model Pembelian Asuransi (Nominal)") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Ringkasan Koefisien Model Pembelian Asuransi (Nominal)
term estimate std.error statistic p.value
(Intercept) -1.773 0.350 -5.063636 0.000
travel_typeInternational 0.775 0.191 4.051494 0.000
risk_levelMedium 0.355 0.223 1.589841 0.112
risk_levelHigh 0.820 0.276 2.973290 0.003
risk_levelVery High 1.819 0.373 4.877162 0.000
trip_duration 0.113 0.036 3.133421 0.002

Interpretasi

Sama dengan interpretasi sebelumnya: variabel perjalanan internasional, risiko tinggi, risiko sangat tinggi, dan durasi semuanya secara signifikan meningkatkan kemungkinan membeli asuransi.

Pemilihan Model Regresi Logistik dan Evaluasi

Membangun Model Regresi Logistik: Pendekatan Confrmatory dan Exploratory

Simulasi Data

Data terdiri dari 300 observasi dengan tiga prediktor: tingkat pendidikan (dalam tahun), status tempat tinggal (kota atau desa), dan tingkat kepercayaan terhadap pemerintah (skala 1–10). Variabel respons adalah kepatuhan memakai masker (1 = patuh, 0 = tidak), yang dibentuk berdasarkan kombinasi logit dari ketiga prediktor. Simulasi ini digunakan untuk mengevaluasi model regresi logistik dan membandingkan kinerjanya dengan berbagai metode seleksi dan evaluasi model.

library(knitr)
library(dplyr)
library(ggplot2)
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
## The following objects are masked from 'package:DescTools':
## 
##     MAE, RMSE
library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
library(DescTools)
set.seed(123)

n <- 300

# x1: Tingkat edukasi (rata-rata sekitar 14 tahun pendidikan formal)
x1 <- rnorm(n, mean = 14, sd = 2)

# x2: Status tinggal di kota besar (1 = kota, 0 = desa)
x2 <- rbinom(n, 1, 0.6)

# x3: Kepercayaan terhadap pemerintah (skala 1–10)
x3 <- round(runif(n, min = 1, max = 10), 1)

# Model logit: konstruksi logika dari pengaruh faktor-faktor terhadap kepatuhan masker
lin_pred <- -2 + 0.25 * x1 + 1 * x2 + 0.3 * x3

# Probabilitas menggunakan masker
p <- 1 / (1 + exp(-lin_pred))

# y: Kepatuhan (1 = patuh, 0 = tidak)
y <- rbinom(n, 1, p)

# Gabungkan ke data frame
df <- data.frame(y = as.factor(y), x1, x2, x3)
head(df)

Pemilihan Model

model_full <- glm(y ~ x1 + x2 + x3, data = df, family = binomial)
summary(model_full)
## 
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  2.88296    2.55022   1.130    0.258
## x1          -0.04047    0.17207  -0.235    0.814
## x2           0.87110    0.66178   1.316    0.188
## x3           0.12069    0.13544   0.891    0.373
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 87.687  on 299  degrees of freedom
## Residual deviance: 85.267  on 296  degrees of freedom
## AIC: 93.267
## 
## Number of Fisher Scoring iterations: 6

Interpretasi

  • Koefisien tidak signifikan (p-value semua > 0.05), artinya tidak ada prediktor yang secara signifikan mempengaruhi kepatuhan masker.
  • AIC = 93.27, digunakan untuk perbandingan model.

Metode Stepwise: Forward, Backward, dan Kedua Arah

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

Interpretasi

  • AIC model stepwise lebih rendah (89.68) dibandingkan model penuh (93.27), artinya model stepwise lebih baik.

Evaluasi Model: ROC dan AUC

pred_prob <- predict(step_both, type = "response")
roc_obj <- roc(df$y, pred_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj, main = "Kurva ROC", col = "blue")

Interpretasi

  • Kurva ROC menggambarkan hubungan antara True Positive Rate (Sensitivity) di sumbu Y dan False Positive Rate (1 - Specificity) di sumbu X untuk berbagai nilai ambang (threshold).

  • Titik (0,1) pada grafik adalah titik ideal (100% sensitif dan 100% spesifik)

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

Interpretasi

  • AUC = 0.5, artinya model tidak memiliki kemampuan klasifikasi yang lebih baik dari tebakan acak.

Pseudo R-Squared

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

Interpretasi

  • Semua nilai = 0 → model sangat buruk dalam menjelaskan variasi variabel respons.

Tabel Klasifkasi dan Evaluasi

pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
conf_matrix <- confusionMatrix(factor(pred_class), df$y, positive = "1")
## Warning in confusionMatrix.default(factor(pred_class), df$y, positive = "1"):
## Levels are not in the same order for reference and data. Refactoring data to
## match.
conf_matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0   0   0
##          1  10 290
##                                           
##                Accuracy : 0.9667          
##                  95% CI : (0.9396, 0.9839)
##     No Information Rate : 0.9667          
##     P-Value [Acc > NIR] : 0.583052        
##                                           
##                   Kappa : 0               
##                                           
##  Mcnemar's Test P-Value : 0.004427        
##                                           
##             Sensitivity : 1.0000          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.9667          
##          Neg Pred Value :    NaN          
##              Prevalence : 0.9667          
##          Detection Rate : 0.9667          
##    Detection Prevalence : 1.0000          
##       Balanced Accuracy : 0.5000          
##                                           
##        'Positive' Class : 1               
##

Interpretasi

  • Model hanya memprediksi kelas 1 (patuh).
  • Semua data kelas 1 (290) diprediksi benar → Sensitivity = 1
  • Semua data kelas 0 (10) diprediksi salah → Specificity = 0
  • Akurasi = 96.67% terlihat tinggi, tapi Kappa = 0 menunjukkan tidak ada kesesuaian selain kebetulan.
  • Balanced Accuracy: 0.5 → Seimbang hanya secara matematis
  • Mcnemar’s Test: p = 0.004 → Perbedaan signifikan antara kesalahan klasifikasi
conf_matrix$byClass[c("Sensitivity", "Specificity")]
## Sensitivity Specificity 
##           1           0

Interpretasi

  • Sensitivitas = 1, Spesifisitas = 0 → model hanya mampu mendeteksi kelas mayoritas.

Metode Perbandingan Model dalam Regresi Logistik

library(MASS)
library(broom)
library(DescTools)

Simulasi Data

Data terdiri dari 300 responden, dengan variabel respons berupa status kepemilikan asuransi (1 = memiliki, 0 = tidak). Variabel prediktor terdiri dari usia (x1), status pekerjaan (x2), dan pendapatan bulanan (x3).

set.seed(456)

n <- 300

# x1: Usia (20–60 tahun)
x1 <- round(runif(n, 20, 60), 0)

# x2: Status pekerjaan (1 = tetap, 0 = tidak tetap)
x2 <- rbinom(n, 1, 0.6)

# x3: Pendapatan bulanan (rata-rata 5 juta, sd 2 juta)
x3 <- round(rnorm(n, mean = 5, sd = 2), 2)
x3[x3 < 0] <- 0.5  # pendapatan tidak boleh negatif

# Linear predictor
lin_pred <- -3 + 0.06 * x1 + 1.2 * x2 + 0.4 * x3

# Probabilitas memiliki asuransi
p <- 1 / (1 + exp(-lin_pred))

# Respon: y = memiliki asuransi (1) atau tidak (0)
y <- rbinom(n, 1, p)

# Data frame
data <- data.frame(y = as.factor(y), x1, x2, x3)
head(data)

Pembuatan Model

# Model 1: hanya usia
model1 <- glm(y ~ x1, data = data, family = binomial)

# Model 2: usia + status pekerjaan
model2 <- glm(y ~ x1 + x2, data = data, family = binomial)

# Model 3: semua variabel
model3 <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)

Perbandingan AIC dan Deviance

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

Interpretasi

  • Penambahan variabel menurunkan AIC dan deviance.
  • Model 3 terbaik dengan AIC = 224.53.

Likelihood-Ratio Test

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

Interpretasi

  • Setiap penambahan variabel memberikan peningkatan signifikan terhadap model (p < 0.05).

Prinsip Parsimony

Evaluasi Tabel Klasifkasi dan Akurasi Model

pred_prob <- predict(model3, type = "response")
pred_class <- factor(ifelse(pred_prob >= 0.5, 1, 0))
conf_matrix <- confusionMatrix(pred_class, data$y, positive = "1")
conf_matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0  13   7
##          1  37 243
##                                           
##                Accuracy : 0.8533          
##                  95% CI : (0.8082, 0.8914)
##     No Information Rate : 0.8333          
##     P-Value [Acc > NIR] : 0.1984          
##                                           
##                   Kappa : 0.3053          
##                                           
##  Mcnemar's Test P-Value : 1.232e-05       
##                                           
##             Sensitivity : 0.9720          
##             Specificity : 0.2600          
##          Pos Pred Value : 0.8679          
##          Neg Pred Value : 0.6500          
##              Prevalence : 0.8333          
##          Detection Rate : 0.8100          
##    Detection Prevalence : 0.9333          
##       Balanced Accuracy : 0.6160          
##                                           
##        'Positive' Class : 1               
##

Interpretasi

  • Akurasi: 85.33% → Cukup tinggi
  • Kappa: 0.305 → Kesesuaian sedang
  • Balanced Accuracy: 0.616 → Kinerja total terhadap kedua kelas cukup baik
  • PPV (Precision): 0.8679 → Prediksi positif sebagian besar benar
  • NPV: 0.65 → Prediksi negatif masih bisa dipercaya, tapi tidak kuat
  • Mcnemar’s p-value: < 0.001 → Ada ketidakseimbangan prediksi FN dan FP
  • Kesimpulan: Model sangat baik dalam mengenali kelas 1, tapi masih lemah dalam deteksi kelas 0.

Sensitivitas dan Spesifsitas

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

Interpretasi

  • Sensitivity: 0.972 → Model sangat baik mendeteksi yang memiliki asuransi

  • Specificity: 0.260 → Lemah dalam mengenali yang tidak memiliki

Detail ROCPenjelasan Kurva ROC (Receiver Operating Characteristic)

Visualisasi dalam R

library(pROC)
set.seed(123)
x1 <- rnorm(200)
x2 <- rbinom(200, 1, 0.5)
x3 <- rnorm(200)
lin_pred <- -1 + 1.5 * x1 - 0.7 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(200, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
pred <- predict(model, type = "response")
roc_obj <- roc(data$y, pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj)

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

Interpretasi

  • AUC = 0.868 → model memiliki performa klasifikasi yang sangat baik.

Simulasi Pemilihan Threshold Optimal

thresholds <- seq(0.1, 0.9, by = 0.05)
results <- data.frame(Threshold = thresholds)
results$Sensitivity <- sapply(thresholds, function(t) {
pred_class <- ifelse(pred >= t, 1, 0)
cm <- table(Pred = pred_class, Obs = data$y)
TP <- cm["1", "1"]
FN <- cm["0", "1"]
TP / (TP + FN)
})
results$Specificity <- sapply(thresholds, function(t) {
pred_class <- ifelse(pred >= t, 1, 0)
cm <- table(Pred = pred_class, Obs = data$y)
TN <- cm["0", "0"]
FP <- cm["1", "0"]
TN / (TN + FP)
})
print(results)
##    Threshold Sensitivity Specificity
## 1       0.10  0.91489362   0.5947712
## 2       0.15  0.85106383   0.6862745
## 3       0.20  0.80851064   0.7320261
## 4       0.25  0.76595745   0.7712418
## 5       0.30  0.72340426   0.8104575
## 6       0.35  0.68085106   0.8366013
## 7       0.40  0.61702128   0.8954248
## 8       0.45  0.59574468   0.9150327
## 9       0.50  0.51063830   0.9281046
## 10      0.55  0.51063830   0.9477124
## 11      0.60  0.42553191   0.9607843
## 12      0.65  0.36170213   0.9738562
## 13      0.70  0.29787234   0.9803922
## 14      0.75  0.19148936   0.9869281
## 15      0.80  0.12765957   0.9869281
## 16      0.85  0.06382979   1.0000000
## 17      0.90  0.02127660   1.0000000

Interpretasi

  • Trade-off antara sensitivitas dan spesifisitas tergantung ambang batas klasifikasi.

Precision-Recall Curve (PR Curve)

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
## The following object is masked from 'package:magrittr':
## 
##     set_names
set.seed(123)
x1 <- rnorm(200)
x2 <- rbinom(200, 1, 0.5)
x3 <- rnorm(200)
lin_pred <- -1 + 1.5 * x1 - 0.7 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(200, 1, p)
data <- data.frame(y = y, x1, x2, x3)
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
prob <- predict(model, type = "response")
pr <- pr.curve(scores.class0 = prob[data$y == 1],
scores.class1 = prob[data$y == 0],
curve = TRUE)
plot(pr)

Interpretasi

  • PR Curve cocok jika data tidak seimbang (class imbalance).

Pseudo R-squared pada Regresi Logistik

Simulasi Data

Data terdiri dari 300 pelanggan dengan tiga prediktor: kecepatan internet (x1), kualitas layanan pelanggan (x2), dan harga langganan bulanan (x3). Variabel respons berupa kepuasan pelanggan (1 = puas, 0 = tidak puas) dibentuk dari kombinasi logit terhadap ketiga prediktor. Simulasi ini digunakan untuk menghitung nilai pseudo R-squared sebagai ukuran kebaikan model regresi logistik, dengan membandingkan model penuh terhadap model null.

set.seed(789)

n <- 300

# x1: Kecepatan internet (10–100 Mbps)
x1 <- round(runif(n, 10, 100), 1)

# x2: Layanan pelanggan memuaskan (1 = ya, 0 = tidak)
x2 <- rbinom(n, 1, 0.7)

# x3: Harga langganan bulanan (200–800 ribu rupiah)
x3 <- round(runif(n, 2, 8), 1)

# Linear predictor
lin_pred <- -4 + 0.05 * x1 + 1.5 * x2 - 0.8 * x3

# Probabilitas puas
p <- 1 / (1 + exp(-lin_pred))

# y: Kepuasan pelanggan (1 = puas, 0 = tidak puas)
y <- rbinom(n, 1, p)

# Data frame
data <- data.frame(y = as.factor(y), x1 = x1, x2 = x2, x3 = x3)

Model Logistik dan Null Model

model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
model_null <- glm(y ~ 1, data = data, family = binomial)

Perhitungan Manual R-squared

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

Perhitungan Otomatis dengan Package Tambahan

Menggunakan pscl

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(model)
## fitting null model for pseudo-r2
##         llh     llhNull          G2    McFadden        r2ML        r2CU 
## -37.1896054 -59.5545730  44.7299353   0.3755374   0.1385169   0.4227115

Menggunakan rcompanion

library(rcompanion)
nagelkerke(model)
## $Models
##                                               
## Model: "glm, y ~ x1 + x2 + x3, binomial, data"
## Null:  "glm, y ~ 1, binomial, data"           
## 
## $Pseudo.R.squared.for.model.vs.null
##                              Pseudo.R.squared
## McFadden                             0.375537
## Cox and Snell (ML)                   0.138517
## Nagelkerke (Cragg and Uhler)         0.422712
## 
## $Likelihood.ratio.test
##  Df.diff LogLik.diff Chisq   p.value
##       -3     -22.365 44.73 1.056e-09
## 
## $Number.of.observations
##           
## Model: 300
## Null:  300
## 
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
## 
## $Warnings
## [1] "None"

Menggunakan DescTools

library(DescTools)
PseudoR2(model, which = "all")
##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.3755374       0.3083721       0.1385169       0.4227115       0.1297536 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.4565636       0.2262188       0.6818220       0.2400787      82.3792107 
##             BIC          logLik         logLik0              G2 
##      97.1943406     -37.1896054     -59.5545730      44.7299353

Interpretasi

Model regresi logistik yang digunakan menunjukkan kinerja yang cukup baik. Nilai McFadden R² sebesar 0,375 dan Nagelkerke R² sebesar 0,423 menunjukkan bahwa model memiliki kemampuan penjelasan yang cukup kuat terhadap variabel dependen. Selain itu, uji likelihood ratio signifikan (p < 0,001), yang berarti model secara signifikan lebih baik dibandingkan model tanpa prediktor (null model).

Apa itu Distribusi Multinomial

Distribusi multinomial adalah perluasan dari distribusi binomial untuk lebih dari dua kategori.

Jika \(X_1, X_2, ..., X_k\) menyatakan banyaknya kejadian dalam masing-masing dari \(k\) kategori, maka:

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

dengan \(\sum_{i=1}^{k} x_i = n\) dan \(\sum_{i=1}^{k} p_i = 1\).

Studi Kasus

Sebuah survei dilakukan terhadap 10 orang tentang genre musik favorit mereka:

  • Pop: 5
  • Rock: 3
  • Jazz: 2

Probabilitas teori:

  • \(p_{Pop} = 0.4\)
  • \(p_{Rock} = 0.4\)
  • \(p_{Jazz} = 0.2\)

Hitung Peluang Manual di R:

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

# Komponen koefisien
faktorial_total <- factorial(n)
faktorial_x <- prod(factorial(x))
koefisien <- faktorial_total / faktorial_x

# Hitung peluang
peluang <- koefisien * prod(p^x)
peluang
## [1] 0.06606029

Multinomial Logistic Regression

Model ini digunakan untuk memodelkan hubungan antara satu variabel respon kategorik (>2 kategori) dan satu atau lebih variabel prediktor.

Misalkan \(Y\) memiliki \(K\) kategori, dan kita pilih referensi (baseline) kategori \(K\), maka model logit untuk kategori \(j\) adalah:

\[ \log\left(\frac{P(Y = j)}{P(Y = K)}\right) = \beta_{j0} + \beta_{j1}x_1 + \cdots + \beta_{jp}x_p, \]

untuk \(j = 1, 2, \dots, K-1\).

Baseline-category logit model

Baseline-category logit model adalah model regresi logistik untuk variabel respon kategorik dengan lebih dari dua kategori (nominal). Model ini menggunakan satu kategori sebagai acuan (baseline) dan membandingkan kategori lainnya terhadap baseline tersebut dalam bentuk logit.

\[ \log\left(\frac{\pi_j}{\pi_c}\right), \quad j = 1, \dots, c-1 \]

dengan:

  • \(\pi_j\) adalah probabilitas respon berada di kategori \(j\)
  • \(\pi_c\) adalah probabilitas respon berada di kategori acuan (baseline)

Terdapat sebanyak \((c - 1)\) fungsi logit.

Catatan: Kategori baseline bisa ditentukan secara eksplisit, tetapi default di R adalah kategori terakhir.

Jika terdapat satu prediktor \(x\), maka bentuk umum model logit-nya adalah:

\[ \log\left(\frac{\pi_j}{\pi_c}\right) = \alpha_j + \beta_j x, \quad j = 1, \dots, c-1 \]

Contoh Kasus: 3 Kategori Respon

Misalkan respon \(Y\) memiliki tiga kategori: \(Y \in \{1, 2, 3\}\), dan kita gunakan kategori ke-3 sebagai baseline. Maka:

\[ \log\left(\frac{\pi_1}{\pi_3}\right) = \alpha_1 + \beta_1 x \]

\[ \log\left(\frac{\pi_2}{\pi_3}\right) = \alpha_2 + \beta_2 x \]

Relasi antar kategori selain baseline:

\[ \log\left(\frac{\pi_1}{\pi_2}\right) = \log\left(\frac{\pi_1/\pi_3}{\pi_2/\pi_3}\right) = \log\left(\frac{\pi_1}{\pi_3}\right) - \log\left(\frac{\pi_2}{\pi_3}\right) = (\alpha_1 - \alpha_2) + (\beta_1 - \beta_2)x \]

Estimasi Parameter

Estimasi dilakukan dengan metode Maximum Likelihood menggunakan algoritma iteratif (seperti Newton-Raphson).

Log-likelihood:

\[ \ell(\beta) = \sum_{i=1}^{n} \sum_{j=1}^{K} y_{ij} \log(\pi_{ij}) \]

dengan: - \(\pi_{ij} = P(Y_i = j | x_i)\) - \(y_{ij} = 1\) jika \(Y_i = j\)

Contoh Kasus

Simulasi Data

library(nnet)
set.seed(123)

n <- 150
Genre <- sample(c("Pop", "Jazz", "Rock"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 30, sd = 5))
ListeningTime <- round(pmax(rnorm(n, mean = 10, sd = 3), 0))

Platform <- sapply(Genre, function(g) {
  if (g == "Pop") {
    sample(c("Spotify", "AppleMusic", "YouTube"), size = 1, prob = c(0.5, 0.3, 0.2))
  } else if (g == "Jazz") {
    sample(c("Spotify", "AppleMusic", "YouTube"), size = 1, prob = c(0.2, 0.5, 0.3))
  } else {
    sample(c("Spotify", "AppleMusic", "YouTube"), size = 1, prob = c(0.3, 0.3, 0.4))
  }
})

df <- data.frame(Platform = factor(Platform), Age, Genre = factor(Genre), ListeningTime)
df$Platform <- relevel(df$Platform, ref = "Spotify") # baseline
head(df)

Estimasi Model

model <- multinom(Platform ~ Age + ListeningTime + Genre, data = df)
## # weights:  18 (10 variable)
## initial  value 164.791843 
## iter  10 value 155.170878
## final  value 155.008122 
## converged
summary(model)
## Call:
## multinom(formula = Platform ~ Age + ListeningTime + Genre, data = df)
## 
## Coefficients:
##            (Intercept)          Age ListeningTime  GenrePop  GenreRock
## AppleMusic   1.6287355 -0.006700859   -0.03947295 -1.601187 -0.9527845
## YouTube      0.5124881  0.041165907   -0.09094170 -2.126068 -0.7515400
## 
## Std. Errors:
##            (Intercept)        Age ListeningTime  GenrePop GenreRock
## AppleMusic    1.381306 0.03900864    0.06746306 0.5346170 0.5181676
## YouTube       1.460503 0.04227322    0.07202683 0.6125704 0.5299225
## 
## Residual Deviance: 310.0162 
## AIC: 330.0162

Nilai P-Value dan Interpretasi

# Koefisien
coef(model)
##            (Intercept)          Age ListeningTime  GenrePop  GenreRock
## AppleMusic   1.6287355 -0.006700859   -0.03947295 -1.601187 -0.9527845
## YouTube      0.5124881  0.041165907   -0.09094170 -2.126068 -0.7515400
# P-value
z <- summary(model)$coefficients / summary(model)$standard.errors
p <- 2 * (1 - pnorm(abs(z)))
p
##            (Intercept)       Age ListeningTime     GenrePop  GenreRock
## AppleMusic   0.2383475 0.8636114     0.5584773 0.0027442927 0.06595088
## YouTube      0.7256645 0.3301530     0.2067298 0.0005190413 0.15613021

Interpretasi

  • Variabel Genre Pop dan Rock berpengaruh signifikan terhadap pemilihan platform AppleMusic dibanding Spotify, karena nilai p-value untuk GenrePop adalah 0.0027 dan GenreRock sebesar 0.0659.
  • Untuk platform YouTube, genre Pop juga signifikan dibanding Spotify (p-value = 0.0005).
  • Namun, variabel Age dan ListeningTime tidak signifikan terhadap pemilihan platform karena p-value > 0.05.

Artinya, preferensi genre lebih menentukan pemilihan platform musik daripada usia atau lama waktu mendengarkan.

Prediksi dan Validasi

df$Predicted <- predict(model)
table(Predicted = df$Predicted, Actual = df$Platform)
##             Actual
## Predicted    Spotify AppleMusic YouTube
##   Spotify         28         14      10
##   AppleMusic      12         29      17
##   YouTube          8         13      19

Interpretasi

Tabel klasifikasi menunjukkan seberapa akurat model dalam memprediksi platform musik. Contoh:

  • Dari pengguna yang menggunakan Spotify, 28 diprediksi benar, sisanya salah klasifikasi.
  • Untuk AppleMusic, 29 diprediksi benar dari total 60.
  • Untuk YouTube, prediksi tepat sebanyak 19 dari 46.

Hal ini menunjukkan bahwa model cukup baik memprediksi, terutama untuk kategori AppleMusic dan YouTube, meskipun masih terdapat mis-klasifikasi yang cukup banyak.

library(ggplot2)

ggplot(df, aes(x = Age, y = ListeningTime, color = Predicted)) +
  geom_point(size = 2) +
  labs(title = "Multinomial Logistic Regression Predictions",
       x = "Age", y = "Listening Time") +
  theme_minimal()

Interpretasi

Grafik menunjukkan hasil prediksi model regresi logistik multinomial terhadap preferensi platform musik berdasarkan usia dan waktu mendengarkan. Secara umum, Apple Music lebih dipilih oleh pengguna berusia muda dengan waktu mendengarkan tinggi, Spotify dominan pada usia menengah, dan YouTube cenderung dipilih oleh pengguna berusia lebih tua. Meskipun pola klasifikasi cukup terlihat, masih terdapat tumpang tindih antar kategori, yang menunjukkan bahwa prediksi belum sepenuhnya akurat hanya dengan dua variabel ini.

Regresi logistik ordinal

Regresi logistik ordinal digunakan ketika variabel respon \(Y\) bersifat ordinal (memiliki urutan), misalnya tingkat kepuasan: Rendah, Sedang, Tinggi. Model ini berbeda dengan: - Regresi logistik biner: hanya 2 kategori - Regresi logistik multinomial: kategori lebih dari 2 tetapi tidak berurutan

Konsep Cumulative Logit Model

Model yang digunakan adalah Cumulative Logit Model dengan asumsi proportional odds:

\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j + \beta x \]

dengan: - \(\alpha_j\): intercept khusus untuk kategori ke-\(j\) - \(\beta\): koefisien regresi (sama untuk semua kategori kumulatif)

Untuk \(c\) kategori, terdapat \((c - 1)\) model logit kumulatif.

Interpretasi Koefisien

Koefisien \(\beta\) menjelaskan efek \(x\) terhadap kemungkinan berada pada kategori yang lebih rendah atau sama.

  • Jika \(\beta > 0\): semakin besar \(x\), semakin tinggi peluang berada di kategori rendah.
  • Jika \(\beta < 0\): semakin besar \(x\), semakin besar peluang berada di kategori tinggi.

Odds ratio: \[ \text{OR} = e^\beta \]

Contoh Data

set.seed(2025)
n <- 250

# Variabel prediktor
tasks_per_week <- round(runif(n, 3, 15))         # Jumlah tugas per minggu
sleep_hours <- round(rnorm(n, mean = 6, sd = 1), 1)  # Jam tidur per malam

# Kombinasi linear (logit ordinal)
linpred <- -0.4 * sleep_hours + 0.3 * tasks_per_week + rnorm(n, 0, 1)

# Konversi menjadi tingkat stres ordinal
stress <- cut(linpred,
              breaks = c(-Inf, 3, 6, Inf),
              labels = c("Rendah", "Sedang", "Tinggi"),
              ordered_result = TRUE)

# Dataset
df_stress <- data.frame(stress, tasks_per_week, sleep_hours)
head(df_stress)

Estimasi Model Ordinal

model_ord <- polr(stress ~ tasks_per_week + sleep_hours, data = df_stress, Hess = TRUE)
(ctable <- coef(summary(model_ord)))
##                     Value Std. Error   t value
## tasks_per_week  0.2809024  0.1060779  2.648078
## sleep_hours    -0.2506101  0.2808765 -0.892243
## Rendah|Sedang   4.4097390  2.0633794  2.137144
## Sedang|Tinggi  84.8715804  2.0633794 41.132319

Nilai P-Value

(ctable <- coef(summary(model_ord)))
##                     Value Std. Error   t value
## tasks_per_week  0.2809024  0.1060779  2.648078
## sleep_hours    -0.2506101  0.2808765 -0.892243
## Rendah|Sedang   4.4097390  2.0633794  2.137144
## Sedang|Tinggi  84.8715804  2.0633794 41.132319
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = round(p, 4)))
##                     Value Std. Error   t value p value
## tasks_per_week  0.2809024  0.1060779  2.648078  0.0081
## sleep_hours    -0.2506101  0.2808765 -0.892243  0.3723
## Rendah|Sedang   4.4097390  2.0633794  2.137144  0.0326
## Sedang|Tinggi  84.8715804  2.0633794 41.132319  0.0000

Prediksi Probabilitas

# Prediksi probabilitas stres untuk kombinasi tugas dan jam tidur
newdata <- data.frame(
  tasks_per_week = c(5, 8, 11),
  sleep_hours = c(7, 6, 5)
)
predict(model_ord, newdata = newdata, type = "probs")
##      Rendah      Sedang Tinggi
## 1 0.9915029 0.008497093      0
## 2 0.9750638 0.024936247      0
## 3 0.9290945 0.070905460      0

Goodness-of-Fit dan Proportional Odds

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

Alternatif Model Ordinal

Selain cumulative logit, terdapat beberapa model ordinal lainnya: - Adjacent-category logit - Continuation-ratio (sequential) logit

Model alternatif digunakan jika asumsi proportional odds tidak terpenuhi.

Kesimpulan

  • Regresi ordinal efektif untuk respon berurutan.
  • Model cumulative logit menginterpretasikan efek dalam bentuk log-odds kumulatif.
  • Implementasi di R dengan fungsi polr() dari package MASS.
  • Untuk validasi lebih lanjut, gunakan uji devians atau likelihood ratio test.

Asumsi Paralelisme dalam Regresi Logistik Ordinal

Model cumulative logit mengasumsikan proportional odds atau parallel lines assumption, yaitu:

\[ \log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j + \beta x, \quad j = 1, \dots, c - 1 \]

  • Hanya \(\alpha_j\) yang berbeda-beda.
  • Koefisien \(\beta\) tetap sama untuk semua fungsi logit kumulatif.

Visualisasi: Kurva logit kumulatif dari tiap kategori terhadap prediktor akan paralel (kemiringan sama, intercept berbeda).

Konsekuensi Pelanggaran Asumsi: - Efek prediktor berbeda untuk tiap batas kategori. - Model cumulative logit tidak valid. - Gunakan model alternatif: Generalized Ordinal Logistic Regression, Partial Proportional Odds Model.

Pengujian Asumsi Paralelisme:

Gunakan Likelihood Ratio Test atau Brant Test.

library(brant)
newdata <- data.frame(
  tasks_per_week = c(5, 8, 11),
  sleep_hours = c(7, 6, 5)
)
predict(model_ord, newdata = newdata, type = "probs")
##      Rendah      Sedang Tinggi
## 1 0.9915029 0.008497093      0
## 2 0.9750638 0.024936247      0
## 3 0.9290945 0.070905460      0

Interpretasi

  • Model cumulative logit menunjukkan bahwa jumlah tugas per minggu adalah prediktor signifikan terhadap tingkat stres mahasiswa.

  • Jam tidur memiliki efek negatif tetapi tidak signifikan secara statistik.

  • Prediksi menunjukkan kecenderungan stres meningkat seiring bertambahnya tugas dan berkurangnya waktu tidur, walaupun seluruh prediksi masih dominan di tingkat “Rendah”.

Log Linear Model

Ringkasan Dalam analisis data kategorik, terdapat beberapa pendekatan statistik yang umum digunakan, antara lain:

  1. Tabel Kontingensi: penyajian frekuensi gabungan dari dua atau lebih variabel kategorik.
  2. Model Loglinear: digunakan untuk memodelkan struktur asosiasi di dalam tabel kontingensi tanpa menganggap ada variabel dependen.
  3. Model Regresi Logistik: digunakan untuk memodelkan probabilitas dari kategori variabel dependen berdasarkan variabel independen.

Meskipun ketiganya dapat digunakan pada data kategorik, pendekatan dan interpretasinya sangat berbeda.

Tabel Kontingensi

Contoh tabel 2x2

tabel_pelatihan <- matrix(c(45, 15, 30, 60), nrow=2, byrow=TRUE,
       dimnames = list(Pelatihan = c("Intensif", "Biasa"),
                       Lulus = c("Ya", "Tidak")))
tabel_pelatihan
##           Lulus
## Pelatihan  Ya Tidak
##   Intensif 45    15
##   Biasa    30    60

Model Loglinear

Model log-linier dua arah (\(I \times J\)) ditulis sebagai:

\[ \log(\mu_{ij}) = \mu + \lambda_i^A + \lambda_j^B + \lambda_{ij}^{AB} \]

library(MASS)
loglm(~ Pelatihan * Lulus, data = tabel_pelatihan)
## Call:
## loglm(formula = ~Pelatihan * Lulus, data = tabel_pelatihan)
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Model Regresi Logistik

Model regresi logistik memodelkan:

\[ \log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 \cdot x \]

data_glm <- data.frame(
  Lulus = c(1, 0, 1, 0),
  Pelatihan = factor(c("Intensif", "Intensif", "Biasa", "Biasa")),
  Frek = c(45, 15, 30, 60)
)

model_logit <- glm(Lulus ~ Pelatihan, weights = Frek, family = binomial, data = data_glm)
summary(model_logit)
## 
## Call:
## glm(formula = Lulus ~ Pelatihan, family = binomial, data = data_glm, 
##     weights = Frek)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -0.6931     0.2236  -3.100  0.00194 ** 
## PelatihanIntensif   1.7918     0.3727   4.808 1.52e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 207.94  on 3  degrees of freedom
## Residual deviance: 182.05  on 2  degrees of freedom
## AIC: 186.05
## 
## Number of Fisher Scoring iterations: 4

Interpretasi

Model Saturated

library(MASS)

# Buat tabel frekuensi (kontingensi)
tabel_pelatihan <- matrix(c(45, 15, 30, 60), nrow=2, byrow=TRUE,
       dimnames = list(Pelatihan = c("Intensif", "Biasa"),
                       Lulus = c("Ya", "Tidak")))

# Fit model saturated (interaksi penuh)
model_saturated <- loglm(~ Pelatihan * Lulus, data = tabel_pelatihan)
summary(model_saturated)
## Formula:
## ~Pelatihan * Lulus
## attr(,"variables")
## list(Pelatihan, Lulus)
## attr(,"factors")
##           Pelatihan Lulus Pelatihan:Lulus
## Pelatihan         1     0               1
## Lulus             0     1               1
## attr(,"term.labels")
## [1] "Pelatihan"       "Lulus"           "Pelatihan:Lulus"
## attr(,"order")
## [1] 1 1 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Interpretasi

Model saturated selalu cocok sempurna dengan data (pasti p = 1), jadi tidak ada deviasi, digunakan sebagai acuan untuk membandingkan model yang lebih sederhana.

Model Independent

Model independen mengasumsikan tidak ada interaksi antar variabel:

\[ \log(\mu_{ij}) = \mu + \lambda_i^A + \lambda_j^B \]

# Contoh data frame
data <- data.frame(
  Pelatihan = c("Intensif", "Intensif", "Biasa", "Biasa"),
  Lulus = c("Ya", "Tidak", "Ya", "Tidak"),
  Frek = c(45, 15, 30, 60)
)

tbl <- xtabs(Frek ~ Pelatihan + Lulus, data = data)

library(MASS)
model_indep <- loglm(~ Pelatihan + Lulus, data = tbl)
summary(model_indep)
## Formula:
## ~Pelatihan + Lulus
## attr(,"variables")
## list(Pelatihan, Lulus)
## attr(,"factors")
##           Pelatihan Lulus
## Pelatihan         1     0
## Lulus             0     1
## attr(,"term.labels")
## [1] "Pelatihan" "Lulus"    
## attr(,"order")
## [1] 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                       X^2 df     P(> X^2)
## Likelihood Ratio 25.89139  1 3.611772e-07
## Pearson          25.00000  1 5.733031e-07

Interpretasi

Model ini tidak cocok dengan data. Nilai p sangat kecil → asumsi independensi ditolak, artinya ada interaksi nyata antara jenis pelatihan dan kelulusan.

Odds Ratio dan Interpretasi

\[ \text{OR} = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}}, \quad \log(\text{OR}) = \log\left(\frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}}\right) \]

  • OR > 1: Pelatihan intensif berhubungan positif dengan kelulusan
  • OR < 1: hubungan negatif
  • OR = 1: tidak ada asosiasi

Estimasi Parameter

tabel <- matrix(c(45, 15, 30, 60), nrow = 2, byrow = TRUE,
                dimnames = list(Pelatihan = c("Intensif", "Biasa"),
                                Lulus = c("Ya", "Tidak")))

logOR <- log((tabel[1,1] * tabel[2,2]) / (tabel[1,2] * tabel[2,1]))
logOR
## [1] 1.791759

Interpretasi

Peserta pelatihan intensif 6 kali lebih besar kemungkinannya untuk lulus dibandingkan peserta pelatihan biasa. Ada asosiasi positif yang kuat.

Model Lebih Sederhana dan Perbandingan Model

anova(model_indep, model_saturated)
## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Pelatihan + Lulus 
## Model 2:
##  ~Pelatihan * Lulus 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   25.89139  1                                    
## Model 2    0.00000  0   25.89139         1              0
## Saturated  0.00000  0    0.00000         0              1

Interpretasi

Model saturated jauh lebih baik dari model independent (ΔX² = 25.89, p = 0), interaksi harus dipertimbangkan.

Studi Kasus: Kepuasan dan Sertifikasi

data_studi <- matrix(c(20, 80,
                       50, 50,
                       40, 60),
                     nrow = 3, byrow = TRUE,
                     dimnames = list(Kepuasan = c("Rendah", "Sedang", "Tinggi"),
                                     Sertifikasi = c("Tidak", "Ya")))
ftable(data_studi)
##          Sertifikasi Tidak Ya
## Kepuasan                     
## Rendah                  20 80
## Sedang                  50 50
## Tinggi                  40 60
loglm(~ Kepuasan + Sertifikasi, data = data_studi)
## Call:
## loglm(formula = ~Kepuasan + Sertifikasi, data = data_studi)
## 
## Statistics:
##                       X^2 df     P(> X^2)
## Likelihood Ratio 20.98240  2 2.777980e-05
## Pearson          20.09569  2 4.327883e-05

Interpretasi

  • Model tanpa interaksi tidak cukup baik, berarti hubungan antara kepuasan dan sertifikasi mungkin bersifat interaktif atau kompleks.

  • Perlu model dengan interaksi untuk menjelaskan data ini dengan baik.

TUGAS MODEL LOG LINEAR 2 ARAH

Model log-linear pada tabel kontingensi.

Model log-linear adalah model yang digunakan untuk menganalisis hubungan antara dua atau lebih variabelkategorik yang disajikan dalam tabel kontingensi. Model ini mengasumsikan bahwa logaritma dari nilai ekspektasi frekuensi sel \((\mu_{ij})\) dapat dinyatakan sebagai penjumlahan efek variabel dan (bila perlu) interaksinya. Untuk tabel 2x2:

\[ \log(\mu_{ij}) = \lambda + \lambda_i^{(A)} + \lambda_j^{(B)} + \lambda_{ij}^{(AB)} \] ## Perbedaan utama antara model log-linear dan model regresi logistik

• Model log-linear digunakan untuk memodelkan frekuensi (count) pada tabel kontingensi dan menguji asosiasi antar variabel kategorik, tanpa menganggap ada variabel respon dan prediktor. • Model regresi logistik digunakan untuk memodelkan probabilitas kejadian suatu outcome (biner) berdasarkan satu atau lebih prediktor (bisa kategorik maupun kontinu).

Analisis Data Tabel Kontingensi 2x2

Misalkan data :

Data:

Sakit Sehat Total
Olahraga: Ya 25 35 60
Olahraga: Tidak 15 25 40
Total 40 60 100

Notasi: - n11 = 25 - n12 = 35 - n21 = 15 - n22 = 25

Model Log Linear

\[ \log(\mu_{ij}) = \lambda + \lambda_i^{(A)} + \lambda_j^{(B)} + \lambda_{ij}^{(AB)} \]

Model tanpa interaksi (independen):

\[ \log(\mu_{ij}) = \lambda + \lambda_i^{(A)} + \lambda_j^{(B)} \]

Estimasi Parameter

Rata-rata log frekuensi sel

\[ \lambda = \frac{1}{4} (\log(25) + \log(35) + \log(15) + \log(25)) = 3.175287 \]

Efek Utama A (Olahraga)

\[ \lambda^{A}_{1} = \frac{1}{2}[(\log(25) + \log(35)) - (\log(15) + \log(25))] / 2 = 0.4236489 \]

\[ \lambda^{A}_{2} = -0.4236489 \]

Efek Utama B (Kesehatan)

\[ \lambda^{B}_{1} = \frac{1}{2}[(\log(25) + \log(35)) - (\log(15) + \log(25))] / 2 = -0.4236489 \] \[ \lambda^{B}_{2} = 0.4236489 \]

Efek Interaksi

\[ \lambda^{AB}_{11} = \frac{1}{4}[\log(25) - \log(35) - \log(15) + \log(25)] = 0.04358835 \]

\[ \lambda^{AB}_{12} = -0.04358835, \quad \lambda^{AB}_{21} = -0.04358835, \quad \lambda^{AB}_{22} = 0.04358835 \]

Ringkasan Parameter

  • \(\lambda = 3.175287\)
  • \(\lambda^A = (0.4236489, - 0.4236489)\)
  • \(\lambda^B = (-0.4236489, 0.4236489)\)
  • \(\lambda^{AB} = \begin{bmatrix} 0.04358835 & -0.04358835 \\ -0.04358835 & 0.04358835 \end{bmatrix}\)

Odds Ratio dan Confidence Interval

Odds Ratio (OR)

\[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} = \frac{25 \times 25}{35 \times 15} = 1.190476 \]

Log Odds Ratio

\[ \log(OR) = \log(1.190476) = 0.1743534 \]

Standard Error (SE)

\[ SE = \sqrt{\frac{1}{25} + \frac{1}{35} + \frac{1}{15} + \frac{1}{25}} = \sqrt{0.2083} = 0.4186145 \]

Interval Kepercayaan 95% untuk log(OR)

\[ 0.1743534 \pm 1.96 \times 0.4186145 = (-0.6461310,0.9948378) \]

Back-transform ke OR

\[ \text{Lower CI} = \exp(0.8968) = 0.5240695, \quad \text{Upper CI} = \exp(2.6868) = 2.7042857 \]

Kesimpulan: Odds Ratio (OR) = 1.190476, dengan 95% CI = (0.5240695, 2.7042857)

Fitting Model Log-Linear dengan R

tabel <- matrix(c(25, 35, 15, 25), nrow=2, byrow=TRUE)
dimnames(tabel) <- list(
  Olahraga = c("Ya", "Tidak"),
  Kesehatan = c("Sakit", "Sehat")
)
tabel
##         Kesehatan
## Olahraga Sakit Sehat
##    Ya       25    35
##    Tidak    15    25
data <- as.data.frame(as.table(tabel))
colnames(data) <- c("Olahraga", "Kesehatan", "Freq")
data
# Model tanpa interaksi
fit_no_inter <- glm(Freq ~ Olahraga + Kesehatan, family = poisson, data = data)
summary(fit_no_inter)
## 
## Call:
## glm(formula = Freq ~ Olahraga + Kesehatan, family = poisson, 
##     data = data)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      3.1781     0.1780  17.859   <2e-16 ***
## OlahragaTidak   -0.4055     0.2041  -1.986    0.047 *  
## KesehatanSehat   0.4055     0.2041   1.986    0.047 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 8.22829  on 3  degrees of freedom
## Residual deviance: 0.17408  on 1  degrees of freedom
## AIC: 26.256
## 
## Number of Fisher Scoring iterations: 3
# Model dengan interaksi
fit_inter <- glm(Freq ~ Olahraga * Kesehatan, family = poisson, data = data)
summary(fit_inter)
## 
## Call:
## glm(formula = Freq ~ Olahraga * Kesehatan, family = poisson, 
##     data = data)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    3.2189     0.2000  16.094   <2e-16 ***
## OlahragaTidak                 -0.5108     0.3266  -1.564    0.118    
## KesehatanSehat                 0.3365     0.2619   1.285    0.199    
## OlahragaTidak:KesehatanSehat   0.1744     0.4186   0.417    0.677    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 8.2283e+00  on 3  degrees of freedom
## Residual deviance: 6.6613e-15  on 0  degrees of freedom
## AIC: 28.082
## 
## Number of Fisher Scoring iterations: 3

Interpretasi

  • Model tanpa interaksi sudah cukup baik.
  • Efek Olahraga dan Kesehatan signifikan secara terpisah, tapi interaksinya tidak signifikan.
  • Artinya, pengaruh olahraga terhadap kesehatan bersifat additif, bukan tergantung satu sama lain.

Analisis Data Tabel Kontingensi 2x3

Sebuah survei dilakukan untuk meneliti hubungan antara Tingkat Pendidikan (SMA, Diploma, Sarjana) dan Status Pekerjaan (Bekerja, Tidak Bekerja). Data dari hasil pengamatan disajikan dalam tabel berikut:

Bekerja Tidak Bekerja
SMA 25 15
Diploma 30 10
Sarjana 28 8

Bentuk umum model log-linear untuk tabel 3x2 (dengan sum-to-zero constraint) adalah sebagai berikut:

\[ \log(\mu_{ij}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_{ij}^{AB} \]

dengan:

  • \(\mu_{ij}\): ekspektasi frekuensi pada baris ke-\(i\), kolom ke-\(j\)
  • A: Tingkat Pendidikan
    (i = 1: SMA, i = 2: Diploma, i = 3: Sarjana)
  • B: Status Pekerjaan
    (j = 1: Bekerja, j = 2: Tidak Bekerja)
  • Constraint: \[ \sum_i \lambda_i^A = 0,\quad \sum_j \lambda_j^B = 0,\quad \sum_i \lambda_{ij}^{AB} = 0,\quad \sum_j \lambda_{ij}^{AB} = 0 \]

Secara eksplisit:

\[ \log(\mu_{ij}) = \lambda + \begin{cases} \lambda_1^A & \text{(SMA)} \\ \lambda_2^A & \text{(Diploma)} \\ \lambda_3^A & \text{(Sarjana)} \end{cases} + \begin{cases} \lambda_1^B & \text{(Bekerja)} \\ \lambda_2^B & \text{(Tidak Bekerja)} \end{cases} + \lambda_{ij}^{AB} \]

Fitting Model Log-Linear di R

# Membuat data frame dari tabel 3x2
tabel3x2 <- matrix(c(25, 15, 30, 10, 28, 8), nrow = 3, byrow = TRUE)
colnames(tabel3x2) <- c("Bekerja", "Tidak Bekerja")
rownames(tabel3x2) <- c("SMA", "Diploma", "Sarjana")
tabel3x2
##         Bekerja Tidak Bekerja
## SMA          25            15
## Diploma      30            10
## Sarjana      28             8
# Ubah ke format data.frame untuk analisis GLM
data3x2 <- as.data.frame(as.table(tabel3x2))
colnames(data3x2) <- c("Pendidikan", "StatusKerja", "Freq")
data3x2
# Model log-linear tanpa interaksi (asumsi independen)
fit_no_inter <- glm(Freq ~ Pendidikan + StatusKerja, family = poisson, data = data3x2)
summary(fit_no_inter)
## 
## Call:
## glm(formula = Freq ~ Pendidikan + StatusKerja, family = poisson, 
##     data = data3x2)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               3.354e+00  1.686e-01  19.893  < 2e-16 ***
## PendidikanDiploma        -4.880e-16  2.236e-01   0.000    1.000    
## PendidikanSarjana        -1.054e-01  2.297e-01  -0.459    0.647    
## StatusKerjaTidak Bekerja -9.223e-01  2.058e-01  -4.482  7.4e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 25.0386  on 5  degrees of freedom
## Residual deviance:  2.4852  on 2  degrees of freedom
## AIC: 38.622
## 
## Number of Fisher Scoring iterations: 4
# Model log-linear dengan interaksi (untuk cek asosiasi)
fit_inter <- glm(Freq ~ Pendidikan * StatusKerja, family = poisson, data = data3x2)
summary(fit_inter)
## 
## Call:
## glm(formula = Freq ~ Pendidikan * StatusKerja, family = poisson, 
##     data = data3x2)
## 
## Coefficients:
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                  3.2189     0.2000  16.094   <2e-16
## PendidikanDiploma                            0.1823     0.2708   0.673    0.501
## PendidikanSarjana                            0.1133     0.2752   0.412    0.680
## StatusKerjaTidak Bekerja                    -0.5108     0.3266  -1.564    0.118
## PendidikanDiploma:StatusKerjaTidak Bekerja  -0.5878     0.4899  -1.200    0.230
## PendidikanSarjana:StatusKerjaTidak Bekerja  -0.7419     0.5171  -1.435    0.151
##                                               
## (Intercept)                                ***
## PendidikanDiploma                             
## PendidikanSarjana                             
## StatusKerjaTidak Bekerja                      
## PendidikanDiploma:StatusKerjaTidak Bekerja    
## PendidikanSarjana:StatusKerjaTidak Bekerja    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.5039e+01  on 5  degrees of freedom
## Residual deviance: 3.7748e-15  on 0  degrees of freedom
## AIC: 40.136
## 
## Number of Fisher Scoring iterations: 3

Interpretasi

  • Model tanpa interaksi (Pendidikan + StatusKerja) sudah cukup baik → tidak perlu interaksi.

  • Tidak terdapat asosiasi signifikan antara tingkat pendidikan dan status pekerjaan berdasarkan data ini.

  • Satu-satunya efek signifikan adalah: orang cenderung lebih banyak bekerja daripada tidak bekerja, terlepas dari tingkat pendidikannya.

Model Log Linear Tiga Arah

Pada pembahasan sebelumnya, kita telah memahami bahwa salah satu tujuan utama dari penyusunan model log-linear adalah untuk mengestimasi parameter-parameter yang menjelaskan hubungan di antara variabel-variabel kategorik.

Pada materi kali ini, kita akan membahas model log-linear yang lebih kompleks, yaitu model log-linear untuk tabel kontingensi tiga arah. Model ini melibatkan tiga variabel kategorik, sehingga kemungkinan interaksi yang dapat terjadi di dalam model pun menjadi lebih banyak. Dalam konteks ini, interaksi paling tinggi yang dapat dimodelkan adalah interaksi tiga arah, yaitu interaksi yang melibatkan ketiga variabel secara bersamaan.

Model Log-Linear untuk Tabel Tiga Arah

Model log-linear yang melibatkan tiga variabel kategorik (misal: X, Y, dan Z) dapat dibangun dalam berbagai bentuk model, tergantung pada tingkat interaksi yang ingin dimasukkan. Berikut adalah beberapa alternatif model log-linear yang umum digunakan:

Model Saturated

Model ini memuat semua kemungkinan interaksi, termasuk interaksi tiga arah (X, Y, dan Z):

\[ \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} \]

2. Model Homogen

Model ini hanya mengakomodasi interaksi dua arah antar variabel tanpa memasukkan interaksi tiga arah:

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

Model Conditional

a. Conditional pada X Memuat interaksi X dengan Y dan X dengan Z:

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

b. Conditional pada Y Memuat interaksi Y dengan X dan Y dengan Z:

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

c. Conditional pada Z Memuat interaksi Z dengan X dan Z dengan Y:

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

Model Joint Independence

a. Independensi antara X dan Y:

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

b. Independensi antara X dan Z:

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

c. Independensi antara Y dan Z:

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

5. Model Tanpa Interaksi

Model ini hanya memasukkan efek utama tanpa interaksi antar variabel:

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

Pengujian Interaksi dalam Model Log-Linear Tiga Arah

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

1. Pengujian Interaksi Tiga Arah (XYZ)

  • Bandingkan model saturated dengan model homogen.
  • Jika perbedaan signifikan → terdapat interaksi tiga arah.

** 2. Pengujian 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.

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

Studi Kasus

Tabel berikut menyajikan data hasil survei nasional yang dilakukan pada tahun 2022 mengenai jenis pekerjaan, jenis kelamin responden, dan preferensi kerja jarak jauh (remote work). Susun dan interpretasikan model log-linear paling sederhana (paling parsimonious) untuk data ini. Jelaskan proses yang Anda lakukan dalam menentukan model terbaik serta asosiasi apa saja yang teridentifikasi. Tunjukkan juga bagaimana nilai yang diprediksi dari model menggambarkan asosiasi tersebut.

Tabel Data Survei

Pekerjaan Jenis Kelamin Setuju Remote Tidak Setuju
Profesional Laki-laki 140 40
Profesional Perempuan 130 50
Administratif Laki-laki 100 70
Administratif Perempuan 120 60
Pekerja Lapangan Laki-laki 60 100
Pekerja Lapangan Perempuan 70 90

Keterangan:

  • Pekerjaan: Profesional, Administratif, Pekerja Lapangan
  • Jenis Kelamin: Laki-laki, Perempuan
  • Preferensi Kerja Jarak Jauh: Setuju (Favor), Tidak Setuju (Oppose)

Analisis Log-Linear untuk Tabel Tiga Arah

Package yang Digunakan

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

Input Data

pekerjaan <- factor(rep(c("Profesional", "Administratif", "Pekerja Lapangan"), each = 4))
jenis_kelamin <- factor(rep(c("Laki-laki", "Perempuan"), each = 2, times = 3))
preferensi <- factor(rep(c("Setuju", "Tidak Setuju"), times = 6))
frekuensi <- c(140, 40, 130, 50, 100, 70, 120, 60, 60, 100, 70, 90)

data <- data.frame(
  Pekerjaan = pekerjaan,
  Jenis_Kelamin = jenis_kelamin,
  Preferensi = preferensi,
  Frekuensi = frekuensi
)

Membentuk Tabel Kontingensi 3 Arah

table3d <- xtabs(Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi, data = data)
ftable(table3d)
##                                Preferensi Setuju Tidak Setuju
## Pekerjaan        Jenis_Kelamin                               
## Administratif    Laki-laki                   100           70
##                  Perempuan                   120           60
## Pekerja Lapangan Laki-laki                    60          100
##                  Perempuan                    70           90
## Profesional      Laki-laki                   140           40
##                  Perempuan                   130           50

UJI MODEL INTERAKSI TIGA ARAH (SATURATED VS HOMOGENOUS)

Penentuan Kategori Referensi

data$Pekerjaan <- relevel(data$Pekerjaan, ref = "Pekerja Lapangan")
data$Jenis_Kelamin <- relevel(data$Jenis_Kelamin, ref = "Perempuan")
data$Preferensi <- relevel(data$Preferensi, ref = "Tidak Setuju")

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

# Model saturated
model_saturated <- glm(Frekuensi ~ Pekerjaan * Jenis_Kelamin * Preferensi,
                       family = poisson(link = "log"),
                       data = data)

summary(model_saturated)
## 
## Call:
## glm(formula = Frekuensi ~ Pekerjaan * Jenis_Kelamin * Preferensi, 
##     family = poisson(link = "log"), data = data)
## 
## Coefficients:
##                                                                Estimate
## (Intercept)                                                     4.49981
## PekerjaanAdministratif                                         -0.40547
## PekerjaanProfesional                                           -0.58779
## Jenis_KelaminLaki-laki                                          0.10536
## PreferensiSetuju                                               -0.25131
## PekerjaanAdministratif:Jenis_KelaminLaki-laki                   0.04879
## PekerjaanProfesional:Jenis_KelaminLaki-laki                    -0.32850
## PekerjaanAdministratif:PreferensiSetuju                         0.94446
## PekerjaanProfesional:PreferensiSetuju                           1.20683
## Jenis_KelaminLaki-laki:PreferensiSetuju                        -0.25951
## PekerjaanAdministratif:Jenis_KelaminLaki-laki:PreferensiSetuju -0.07696
## PekerjaanProfesional:Jenis_KelaminLaki-laki:PreferensiSetuju    0.55676
##                                                                Std. Error
## (Intercept)                                                       0.10541
## PekerjaanAdministratif                                            0.16667
## PekerjaanProfesional                                              0.17638
## Jenis_KelaminLaki-laki                                            0.14530
## PreferensiSetuju                                                  0.15936
## PekerjaanAdministratif:Jenis_KelaminLaki-laki                     0.22817
## PekerjaanProfesional:Jenis_KelaminLaki-laki                       0.25712
## PekerjaanAdministratif:PreferensiSetuju                           0.22449
## PekerjaanProfesional:PreferensiSetuju                             0.23041
## Jenis_KelaminLaki-laki:PreferensiSetuju                           0.22817
## PekerjaanAdministratif:Jenis_KelaminLaki-laki:PreferensiSetuju    0.31835
## PekerjaanProfesional:Jenis_KelaminLaki-laki:PreferensiSetuju      0.33451
##                                                                z value Pr(>|z|)
## (Intercept)                                                     42.689  < 2e-16
## PekerjaanAdministratif                                          -2.433 0.014983
## PekerjaanProfesional                                            -3.332 0.000861
## Jenis_KelaminLaki-laki                                           0.725 0.468366
## PreferensiSetuju                                                -1.577 0.114799
## PekerjaanAdministratif:Jenis_KelaminLaki-laki                    0.214 0.830681
## PekerjaanProfesional:Jenis_KelaminLaki-laki                     -1.278 0.201382
## PekerjaanAdministratif:PreferensiSetuju                          4.207 2.59e-05
## PekerjaanProfesional:PreferensiSetuju                            5.238 1.63e-07
## Jenis_KelaminLaki-laki:PreferensiSetuju                         -1.137 0.255397
## PekerjaanAdministratif:Jenis_KelaminLaki-laki:PreferensiSetuju  -0.242 0.808976
## PekerjaanProfesional:Jenis_KelaminLaki-laki:PreferensiSetuju     1.664 0.096033
##                                                                   
## (Intercept)                                                    ***
## PekerjaanAdministratif                                         *  
## PekerjaanProfesional                                           ***
## Jenis_KelaminLaki-laki                                            
## PreferensiSetuju                                                  
## PekerjaanAdministratif:Jenis_KelaminLaki-laki                     
## PekerjaanProfesional:Jenis_KelaminLaki-laki                       
## PekerjaanAdministratif:PreferensiSetuju                        ***
## PekerjaanProfesional:PreferensiSetuju                          ***
## Jenis_KelaminLaki-laki:PreferensiSetuju                           
## PekerjaanAdministratif:Jenis_KelaminLaki-laki:PreferensiSetuju    
## PekerjaanProfesional:Jenis_KelaminLaki-laki:PreferensiSetuju   .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  1.3683e+02  on 11  degrees of freedom
## Residual deviance: -3.1086e-14  on  0  degrees of freedom
## AIC: 98.675
## 
## Number of Fisher Scoring iterations: 3
exp(coef(model_saturated))
##                                                    (Intercept) 
##                                                     90.0000000 
##                                         PekerjaanAdministratif 
##                                                      0.6666667 
##                                           PekerjaanProfesional 
##                                                      0.5555556 
##                                         Jenis_KelaminLaki-laki 
##                                                      1.1111111 
##                                               PreferensiSetuju 
##                                                      0.7777778 
##                  PekerjaanAdministratif:Jenis_KelaminLaki-laki 
##                                                      1.0500000 
##                    PekerjaanProfesional:Jenis_KelaminLaki-laki 
##                                                      0.7200000 
##                        PekerjaanAdministratif:PreferensiSetuju 
##                                                      2.5714286 
##                          PekerjaanProfesional:PreferensiSetuju 
##                                                      3.3428571 
##                        Jenis_KelaminLaki-laki:PreferensiSetuju 
##                                                      0.7714286 
## PekerjaanAdministratif:Jenis_KelaminLaki-laki:PreferensiSetuju 
##                                                      0.9259259 
##   PekerjaanProfesional:Jenis_KelaminLaki-laki:PreferensiSetuju 
##                                                      1.7450142

Hasil Estimasi Koefisien

Hasil Estimasi Koefisien

Berdasarkan output model saturated yang telah diperoleh:

Term Estimate Std. Error z value Pr(> z
(Intercept) 4.49981 0.10541 42.689 < 2e-16 90.00
PekerjaanAdministratif -0.40547 0.16667 -2.433 0.015 0.67
PekerjaanProfesional -0.58779 0.17638 -3.332 0.0009 0.56
Jenis_KelaminLaki-laki 0.10536 0.14530 0.725 0.468 1.11
PreferensiSetuju -0.25131 0.15936 -1.577 0.115 0.78
PekerjaanAdministratif:Jenis_KelaminLaki-laki 0.04879 0.22817 0.214 0.831 1.05
PekerjaanProfesional:Jenis_KelaminLaki-laki -0.32850 0.25712 -1.278 0.201 0.72
PekerjaanAdministratif:PreferensiSetuju 0.94446 0.22449 4.207 2.59e-05 2.57
PekerjaanProfesional:PreferensiSetuju 1.20683 0.23041 5.238 1.63e-07 3.34
Jenis_KelaminLaki-laki:PreferensiSetuju -0.25951 0.22817 -1.137 0.255 0.77
PekerjaanAdministratif:Jenis_KelaminLaki-laki:PreferensiSetuju -0.07696 0.31835 -0.242 0.809 0.93
PekerjaanProfesional:Jenis_KelaminLaki-laki:PreferensiSetuju 0.55676 0.33451 1.664 0.096 1.75

Interpretasi Koefisien

  • Intercept : Frekuensi dasar responden pada kategori referensi (Pekerja Lapangan, Perempuan, Tidak Setuju) sekitar 90 orang.
  • Pekerjaan : Responden Administratif dan Profesional lebih sedikit secara umum dibanding Pekerja Lapangan, tapi lebih cenderung setuju kerja jarak jauh (interaksi signifikan).
  • Jenis Kelamin : Tidak ada perbedaan signifikan antara laki-laki dan perempuan dalam jumlah responden.
  • Preferensi : Tidak ada perbedaan signifikan frekuensi antara setuju dan tidak setuju pada kategori dasar.
  • Interaksi Pekerjaan-Preferensi : Responden Administratif dan Profesional lebih sering setuju kerja remote dibanding Pekerja Lapangan.
  • Interaksi Jenis Kelamin-Preferensi : Tidak signifikan, jenis kelamin tidak berpengaruh pada preferensi.
  • Interaksi Tiga Arah : Tidak signifikan, model dua arah sudah cukup menjelaskan hubungan antar variabel.

Goodness-of-Fit

  • Nilai deviance residual model saturated adalah 0, menunjukkan model fit dengan sempurna data.

  • ANOVA membandingkan model saturated dengan model homogeneous menunjukkan perbedaan yang tidak signifikan (p > 0.05), mendukung bahwa model dua arah sudah cukup.

Kesimpulan

  • Asosiasi utama dalam data ini adalah hubungan antara jenis pekerjaan dan preferensi kerja remote.
  • Jenis kelamin tidak berperan signifikan dalam preferensi kerja jarak jauh.
  • Model log-linear dua arah memberikan gambaran parsimonious yang baik tanpa kehilangan informasi penting.
  • Prediksi frekuensi dari model menunjukkan bahwa pekerja profesional dan administratif lebih cenderung setuju kerja jarak jauh dibanding pekerja lapangan.

Model Homogenous

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

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

model_homogenous <- glm(Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi +
                        Pekerjaan:Jenis_Kelamin + Pekerjaan:Preferensi + Jenis_Kelamin:Preferensi,
                        family = poisson(link = "log"),
                        data = data)
summary(model_homogenous)
## 
## Call:
## glm(formula = Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi + 
##     Pekerjaan:Jenis_Kelamin + Pekerjaan:Preferensi + Jenis_Kelamin:Preferensi, 
##     family = poisson(link = "log"), data = data)
## 
## Coefficients:
##                                               Estimate Std. Error z value
## (Intercept)                                    4.52896    0.09650  46.932
## PekerjaanAdministratif                        -0.36411    0.13916  -2.616
## PekerjaanProfesional                          -0.76879    0.15263  -5.037
## Jenis_KelaminLaki-laki                         0.04924    0.12421   0.396
## PreferensiSetuju                              -0.31923    0.13149  -2.428
## PekerjaanAdministratif:Jenis_KelaminLaki-laki -0.03024    0.15757  -0.192
## PekerjaanProfesional:Jenis_KelaminLaki-laki    0.04169    0.16039   0.260
## PekerjaanAdministratif:PreferensiSetuju        0.90467    0.15879   5.697
## PekerjaanProfesional:PreferensiSetuju          1.47936    0.16673   8.873
## Jenis_KelaminLaki-laki:PreferensiSetuju       -0.12122    0.13300  -0.911
##                                               Pr(>|z|)    
## (Intercept)                                    < 2e-16 ***
## PekerjaanAdministratif                         0.00888 ** 
## PekerjaanProfesional                          4.73e-07 ***
## Jenis_KelaminLaki-laki                         0.69181    
## PreferensiSetuju                               0.01519 *  
## PekerjaanAdministratif:Jenis_KelaminLaki-laki  0.84780    
## PekerjaanProfesional:Jenis_KelaminLaki-laki    0.79493    
## PekerjaanAdministratif:PreferensiSetuju       1.22e-08 ***
## PekerjaanProfesional:PreferensiSetuju          < 2e-16 ***
## Jenis_KelaminLaki-laki:PreferensiSetuju        0.36208    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 136.8254  on 11  degrees of freedom
## Residual deviance:   4.2543  on  2  degrees of freedom
## AIC: 98.929
## 
## Number of Fisher Scoring iterations: 4

Uji Hipotesis

Hipotesis

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

Hitung Selisih Deviance

Deviance.model <- model_homogenous$deviance - model_saturated$deviance
Deviance.model
## [1] 4.254279

Hitung Derajat Bebas

derajat.bebas <- model_homogenous$df.residual - model_saturated$df.residual
derajat.bebas
## [1] 2

Hitung Chi-Square Tabel

chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel
## [1] 5.991465

Keputusan Uji

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

Rangkuman

  • Hipotesis
    • \(H_0: \lambda^{XYZ} = 0\) (Tidak ada interaksi tiga arah; model yang terbentuk adalah model homogenous)
    • \(H_1: \lambda^{XYZ} \ne 0\) (Ada interaksi tiga arah; model yang terbentuk adalah model saturated)
  • Tingkat Signifikansi
    • \(\alpha = 5\%\)
  • Statistik Uji
    • \(\Delta \text{Deviance} = \text{Deviance model homogenous} - \text{Deviance model saturated}\)
    • \(\Delta \text{Deviance} = 4.254 - 0 = 4.254\)
    • \(df = df_{homogenous} - df_{saturated} = 2 - 0 = 2\)
  • Kriteria Uji
    • Tolak \(H_0\) jika \(\Delta \text{Deviance} > \chi^2_{(0.05, 2)} = 5.991\)
  • Keputusan
    • Karena \(4.254 < 5.991\), maka tidak tolak \(H_0\)
  • Interpretasi
    Pada taraf nyata 5%, belum cukup bukti untuk menolak \(H_0\) sehingga tidak ada interaksi tiga arah antara jenis kelamin, pekerjaan, dan preferensi kerja jarak jauh.

Artinya, model homogenous (yang hanya memasukkan interaksi dua arah) sudah cukup menjelaskan hubungan dalam data.

UJI MODEL INTERAKSI DUAARAH(HOMOGENOUS VS CONDITIONAL ON X)

Model Conditional on X

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

model_conditional_X <- glm(Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi +
                             Pekerjaan:Jenis_Kelamin + Pekerjaan:Preferensi,
                           family = poisson(link = "log"),
                           data = data)
summary(model_conditional_X)
## 
## Call:
## glm(formula = Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi + 
##     Pekerjaan:Jenis_Kelamin + Pekerjaan:Preferensi, family = poisson(link = "log"), 
##     data = data)
## 
## Coefficients:
##                                                 Estimate Std. Error z value
## (Intercept)                                    4.554e+00  9.159e-02  49.722
## PekerjaanAdministratif                        -3.513e-01  1.370e-01  -2.564
## PekerjaanProfesional                          -7.472e-01  1.493e-01  -5.006
## Jenis_KelaminLaki-laki                        -2.805e-16  1.118e-01   0.000
## PreferensiSetuju                              -3.795e-01  1.138e-01  -3.334
## PekerjaanAdministratif:Jenis_KelaminLaki-laki -5.716e-02  1.547e-01  -0.369
## PekerjaanProfesional:Jenis_KelaminLaki-laki   -2.058e-15  1.537e-01   0.000
## PekerjaanAdministratif:PreferensiSetuju        9.056e-01  1.587e-01   5.705
## PekerjaanProfesional:PreferensiSetuju          1.478e+00  1.666e-01   8.870
##                                               Pr(>|z|)    
## (Intercept)                                    < 2e-16 ***
## PekerjaanAdministratif                        0.010356 *  
## PekerjaanProfesional                          5.55e-07 ***
## Jenis_KelaminLaki-laki                        1.000000    
## PreferensiSetuju                              0.000856 ***
## PekerjaanAdministratif:Jenis_KelaminLaki-laki 0.711804    
## PekerjaanProfesional:Jenis_KelaminLaki-laki   1.000000    
## PekerjaanAdministratif:PreferensiSetuju       1.16e-08 ***
## PekerjaanProfesional:PreferensiSetuju          < 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: 136.8254  on 11  degrees of freedom
## Residual deviance:   5.0853  on  3  degrees of freedom
## AIC: 97.76
## 
## Number of Fisher Scoring iterations: 4

Pengujian Selisih Deviance (Conditional on X vs Homogenous)

# Deviance of Model
Deviance.model <- model_conditional_X$deviance - model_homogenous$deviance
Deviance.model
## [1] 0.8310679
# Chi Square tabel dengan alpha = 0.05
derajat.bebas <- 4
derajat.bebas
## [1] 4
chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel
## [1] 9.487729
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0 (model sederhana cukup)", "Tolak H0 (perlu model interaksi)")
Keputusan
## [1] "Terima H0 (model sederhana cukup)"

Interpretasi

Karena nilai Deviance.model sebesar 0.83 lebih kecil dari nilai kritis chi-square tabel sebesar 9.49 (dengan derajat bebas 4 dan tingkat signifikansi 0.05), maka keputusan uji adalah “Terima H0 (model sederhana cukup)”.

Artinya, pada taraf nyata 5%, belum ada bukti yang cukup untuk menolak hipotesis nol. Dengan kata lain, model yang hanya memuat efek utama tanpa interaksi sudah cukup untuk menjelaskan data.

Interaksi antara jenis kelamin dan preferensi, baik secara langsung maupun melalui pekerjaan, tidak signifikan secara statistik. Jadi, hubungan antara jenis kelamin dan preferensi tidak bergantung pada pekerjaan, dan model tanpa interaksi adalah model yang parsimonious dan memadai.

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

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

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

# Conditional Association on Y
model_conditional_Y <- glm(Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi +
                             Pekerjaan:Jenis_Kelamin + Jenis_Kelamin:Preferensi,
                           family = poisson(link = "log"),
                           data = data)

summary(model_conditional_Y)
## 
## Call:
## glm(formula = Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi + 
##     Pekerjaan:Jenis_Kelamin + Jenis_Kelamin:Preferensi, family = poisson(link = "log"), 
##     data = data)
## 
## Coefficients:
##                                                 Estimate Std. Error z value
## (Intercept)                                    4.120e+00  9.658e-02  42.657
## PekerjaanAdministratif                         1.178e-01  1.087e-01   1.084
## PekerjaanProfesional                           1.178e-01  1.087e-01   1.084
## Jenis_KelaminLaki-laki                         6.821e-02  1.356e-01   0.503
## PreferensiSetuju                               4.700e-01  9.014e-02   5.214
## PekerjaanAdministratif:Jenis_KelaminLaki-laki -5.716e-02  1.547e-01  -0.369
## PekerjaanProfesional:Jenis_KelaminLaki-laki    1.376e-11  1.537e-01   0.000
## Jenis_KelaminLaki-laki:PreferensiSetuju       -1.133e-01  1.274e-01  -0.890
##                                               Pr(>|z|)    
## (Intercept)                                    < 2e-16 ***
## PekerjaanAdministratif                           0.278    
## PekerjaanProfesional                             0.278    
## Jenis_KelaminLaki-laki                           0.615    
## PreferensiSetuju                              1.85e-07 ***
## PekerjaanAdministratif:Jenis_KelaminLaki-laki    0.712    
## PekerjaanProfesional:Jenis_KelaminLaki-laki      1.000    
## Jenis_KelaminLaki-laki:PreferensiSetuju          0.374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 136.825  on 11  degrees of freedom
## Residual deviance:  90.082  on  4  degrees of freedom
## AIC: 180.76
## 
## Number of Fisher Scoring iterations: 4
# Deviance of Model
Deviance.model <- model_conditional_Y$deviance - model_homogenous$deviance
Deviance.model
## [1] 85.82817
# Derajat bebas:  
# Interaksi Pekerjaan:Jenis_Kelamin = (3-1)*(2-1) = 2
# Interaksi Jenis_Kelamin:Preferensi = (2-1)*(2-1) = 1
# Total df = 2 + 1 = 3
derajat.bebas <- 3
derajat.bebas
## [1] 3
chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 7.814728
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan
## [1] "Tolak"

interpretasi

Karena nilai Deviance.model sebesar 85.83 jauh lebih besar dari nilai kritis chi-square tabel sebesar 7.81 (dengan derajat bebas 3 dan tingkat signifikansi 0.05), maka keputusan uji adalah “Tolak H0”.

Artinya, pada taraf nyata 5%, terdapat bukti yang cukup untuk menolak hipotesis nol. Dengan kata lain, model dengan interaksi antara pekerjaan dan jenis kelamin serta interaksi antara jenis kelamin dan preferensi secara signifikan lebih baik daripada model tanpa interaksi.

Ini menunjukkan bahwa ada interaksi yang signifikan antara faktor-faktor tersebut, sehingga model yang hanya mengandung efek utama tidak cukup untuk menjelaskan data secara tepat.

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

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

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

Model Conditional on Z

# Conditional Association on Z
model_conditional_Z <- glm(Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi +
                             Pekerjaan:Preferensi + Jenis_Kelamin:Preferensi,
                           family = poisson(link = "log"),
                           data = data)

summary(model_conditional_Z)
## 
## Call:
## glm(formula = Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi + 
##     Pekerjaan:Preferensi + Jenis_Kelamin:Preferensi, family = poisson(link = "log"), 
##     data = data)
## 
## Coefficients:
##                                         Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              4.52918    0.08845  51.204  < 2e-16
## PekerjaanAdministratif                  -0.37949    0.11382  -3.334 0.000856
## PekerjaanProfesional                    -0.74721    0.12796  -5.839 5.24e-09
## Jenis_KelaminLaki-laki                   0.04879    0.09880   0.494 0.621436
## PreferensiSetuju                        -0.32305    0.13049  -2.476 0.013302
## PekerjaanAdministratif:PreferensiSetuju  0.90558    0.15872   5.705 1.16e-08
## PekerjaanProfesional:PreferensiSetuju    1.47810    0.16664   8.870  < 2e-16
## Jenis_KelaminLaki-laki:PreferensiSetuju -0.11333    0.12736  -0.890 0.373553
##                                            
## (Intercept)                             ***
## PekerjaanAdministratif                  ***
## PekerjaanProfesional                    ***
## Jenis_KelaminLaki-laki                     
## PreferensiSetuju                        *  
## PekerjaanAdministratif:PreferensiSetuju ***
## PekerjaanProfesional:PreferensiSetuju   ***
## Jenis_KelaminLaki-laki:PreferensiSetuju    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 136.8254  on 11  degrees of freedom
## Residual deviance:   4.4819  on  4  degrees of freedom
## AIC: 95.157
## 
## Number of Fisher Scoring iterations: 4
# Deviance of Model
Deviance.model <- model_conditional_Z$deviance - model_homogenous$deviance
Deviance.model
## [1] 0.2276212
# Derajat bebas:
# Interaksi Pekerjaan:Preferensi = (3-1)*(2-1) = 2
# Interaksi Jenis_Kelamin:Preferensi = (2-1)*(2-1) = 1
# Total df = 2 + 1 = 3
derajat.bebas <- 3
derajat.bebas
## [1] 3
# Chi-square tabel alpha = 0.05
chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel
## [1] 7.814728
# Keputusan uji
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima H0 (model sederhana cukup)", "Tolak H0 (perlu model interaksi)")
Keputusan
## [1] "Terima H0 (model sederhana cukup)"

Interpretasi

Karena nilai Deviance.model sebesar 0.23 jauh lebih kecil dari nilai kritis chi-square tabel sebesar 7.81 (dengan derajat bebas 3 dan tingkat signifikansi 0.05), maka keputusan uji adalah “Terima H0 (model sederhana cukup)”.

Artinya, pada taraf nyata 5%, belum ada bukti yang cukup untuk menolak hipotesis nol. Dengan kata lain, model yang hanya memuat efek utama tanpa interaksi sudah cukup untuk menjelaskan data.

Interaksi antara pekerjaan dan preferensi serta interaksi antara jenis kelamin dan preferensi tidak signifikan secara statistik. Sehingga, model sederhana tanpa interaksi sudah memadai untuk data ini.

PEMILIHAN MODEL TERBAIK

library(DiagrammeR)

grViz("
digraph model_selection {
  
  graph [layout = dot, rankdir = TB]
  
  node [shape = box, style = filled, fillcolor = lightblue, fontname = Helvetica, fontsize = 12]

  Saturated     [label = 'Model Saturated\\n(Pekerjaan * Jenis Kelamin * Preferensi)']
  Homogenous    [label = 'Model Homogenous\\n(Interaksi Dua Arah)', fillcolor = lightgreen]
  Conditional_X [label = 'Model Conditional on X\\n(Pekerjaan)', fillcolor = lightyellow]
  Conditional_Y [label = 'Model Conditional on Y\\n(Jenis Kelamin)', fillcolor = lightyellow]
  Conditional_Z [label = 'Model Conditional on Z\\n(Preferensi)', fillcolor = lightyellow]
  FinalModel    [label = 'Model Terpilih:\\nConditional on Z\\n(Parsimonious)', fillcolor = orange]
  
  Saturated -> Homogenous [label = 'Tidak signifikan\\n(p > 0.05)']
  Homogenous -> Conditional_X [label = 'p > 0.05']
  Homogenous -> Conditional_Y [label = 'p < 0.05']
  Homogenous -> Conditional_Z [label = 'p > 0.05']
  Conditional_Z -> FinalModel [label = 'Model paling sederhana\\ndan masih fit']
}
")

Model Terbaik

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

# Model Terbaik
bestmodel <- glm(Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi +
                             Pekerjaan:Preferensi + Jenis_Kelamin:Preferensi,
                           family = poisson(link = "log"),
                           data = data)

summary(bestmodel)
## 
## Call:
## glm(formula = Frekuensi ~ Pekerjaan + Jenis_Kelamin + Preferensi + 
##     Pekerjaan:Preferensi + Jenis_Kelamin:Preferensi, family = poisson(link = "log"), 
##     data = data)
## 
## Coefficients:
##                                         Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              4.52918    0.08845  51.204  < 2e-16
## PekerjaanAdministratif                  -0.37949    0.11382  -3.334 0.000856
## PekerjaanProfesional                    -0.74721    0.12796  -5.839 5.24e-09
## Jenis_KelaminLaki-laki                   0.04879    0.09880   0.494 0.621436
## PreferensiSetuju                        -0.32305    0.13049  -2.476 0.013302
## PekerjaanAdministratif:PreferensiSetuju  0.90558    0.15872   5.705 1.16e-08
## PekerjaanProfesional:PreferensiSetuju    1.47810    0.16664   8.870  < 2e-16
## Jenis_KelaminLaki-laki:PreferensiSetuju -0.11333    0.12736  -0.890 0.373553
##                                            
## (Intercept)                             ***
## PekerjaanAdministratif                  ***
## PekerjaanProfesional                    ***
## Jenis_KelaminLaki-laki                     
## PreferensiSetuju                        *  
## PekerjaanAdministratif:PreferensiSetuju ***
## PekerjaanProfesional:PreferensiSetuju   ***
## Jenis_KelaminLaki-laki:PreferensiSetuju    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 136.8254  on 11  degrees of freedom
## Residual deviance:   4.4819  on  4  degrees of freedom
## AIC: 95.157
## 
## Number of Fisher Scoring iterations: 4

INTERPRETASI KOEFISIEN MODEL TERBAIK

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

Interpretasi Koefisien Model Terbaik

  • \(\exp(\beta_{\text{PekerjaanAdministratif}}) = \exp(-0{,}379) = 0{,}684 \rightarrow\) nilai odds
    Tanpa memperhatikan jenis kelamin dan preferensi, peluang individu dengan pekerjaan administratif adalah 0{,}684 kali dibandingkan individu dengan pekerjaan lain (baseline: tidak disebutkan, mungkin kategori “lainnya”).

  • \(\exp(\beta_{\text{PekerjaanProfesional}}) = \exp(-0{,}747) = 0{,}474 \rightarrow\) nilai odds
    Tanpa memperhatikan jenis kelamin dan preferensi, peluang individu dengan pekerjaan profesional adalah 0{,}474 kali dibandingkan kelompok referensi pekerjaan lainnya.

  • \(\exp(\beta_{\text{Jenis\_KelaminLaki-laki}}) = \exp(0{,}049) = 1{,}050 \rightarrow\) nilai odds
    Tanpa memperhatikan pekerjaan dan preferensi, peluang laki-laki adalah 1{,}050 kali dibandingkan perempuan.

  • \(\exp(\beta_{\text{PreferensiSetuju}}) = \exp(-0{,}323) = 0{,}724 \rightarrow\) nilai odds
    Tanpa memperhatikan pekerjaan dan jenis kelamin, peluang individu yang setuju (terhadap sesuatu) adalah 0{,}724 kali dibandingkan individu yang tidak setuju.

  • \(\exp(\beta_{\text{PekerjaanAdministratif:Setuju}}) = \exp(0{,}906) = 2{,}474 \rightarrow\) nilai odds ratio
    Pada kelompok yang setuju, odds individu dengan pekerjaan administratif adalah 2{,}47 kali dibandingkan odds yang sama pada kelompok yang tidak setuju.

  • \(\exp(\beta_{\text{PekerjaanProfesional:Setuju}}) = \exp(1{,}478) = 4{,}385 \rightarrow\) nilai odds ratio
    Pada kelompok yang setuju, odds individu dengan pekerjaan profesional adalah 4{,}385 kali dibandingkan kelompok referensi (misalnya pekerjaan lain) dalam kelompok tidak setuju.

  • \(\exp(\beta_{\text{Laki-laki:Setuju}}) = \exp(-0{,}113) = 0{,}893 \rightarrow\) nilai odds ratio
    Pada kelompok yang setuju, odds laki-laki adalah 0{,}893 kali dibandingkan odds perempuan pada kelompok referensi preferensi.

  • \(\exp(\text{Intercept}) = \exp(4{,}530) = 92{,}683 \rightarrow\) baseline odds
    Ini adalah odds dasar (intercept), yaitu odds untuk kategori referensi dari semua variabel (misalnya: perempuan, pekerjaan referensi, preferensi tidak setuju).

NILAI DUGAAN MODEL TERBAIK

 # Fitted values dari model terbaik
 data.frame(
  Pekerjaan = data$Pekerjaan,
  Jenis_Kelamin = data$Jenis_Kelamin,
  Preferensi = data$Preferensi,
  Frekuensi = data$Frekuensi,
  Fitted = bestmodel$fitted.values
)

Contoh Kasus : Analisis Prediktif Diabetes Menggunakan Analisis Regresi Logistik : Pemilihan Model dan Evaluasi Kinerja

Pendahuluan

Diabetes mellitus merupakan salah satu penyakit kronis yang tingkat prevalensinya semakin meningkat secara global. Deteksi dini terhadap risiko diabetes sangat penting untuk mencegah komplikasi jangka panjang. Dalam konteks ini, metode statistik seperti regresi logistik digunakan untuk membangun model prediktif yang dapat membantu mengidentifikasi individu dengan risiko tinggi terhadap diabetes berdasarkan indikator medis.

Urgensi

Penggunaan model prediktif dalam bidang medis dapat meningkatkan efisiensi diagnosis, membantu pengambilan keputusan klinis, dan menekan biaya pemeriksaan. Pemodelan yang akurat dan terukur, serta evaluasi performa model dengan teknik seperti AUC, precision-recall, dan pseudo R-squared menjadi kunci dalam menjamin efektivitas sistem prediksi ini.

Dataset

Dataset yang digunakan adalah “diabetes.csv” dari lokal komputer, memuat 768 observasi dan beberapa variabel prediktor: - Pregnancies: jumlah kehamilan - Glucose: kadar glukosa plasma - BloodPressure: tekanan darah diastolik - SkinThickness: ketebalan lipatan kulit - Insulin: kadar insulin serum - BMI: indeks massa tubuh - DiabetesPedigreeFunction: riwayat keluarga diabetes - Age: usia - Outcome: status diabetes (1 = diabetes, 0 = tidak diabetes)

Analisis dan Evaluasi Model

data <- read.csv("C:/Users/asus/Downloads/diabetes.csv")
data$Outcome <- as.factor(data$Outcome)
head(data)

Model Penuh

model_full <- glm(Outcome ~ ., data = data, family = binomial)
summary(model_full)
## 
## Call:
## glm(formula = Outcome ~ ., family = binomial, data = data)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -8.4046964  0.7166359 -11.728  < 2e-16 ***
## Pregnancies               0.1231823  0.0320776   3.840 0.000123 ***
## Glucose                   0.0351637  0.0037087   9.481  < 2e-16 ***
## BloodPressure            -0.0132955  0.0052336  -2.540 0.011072 *  
## SkinThickness             0.0006190  0.0068994   0.090 0.928515    
## Insulin                  -0.0011917  0.0009012  -1.322 0.186065    
## BMI                       0.0897010  0.0150876   5.945 2.76e-09 ***
## DiabetesPedigreeFunction  0.9451797  0.2991475   3.160 0.001580 ** 
## Age                       0.0148690  0.0093348   1.593 0.111192    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 993.48  on 767  degrees of freedom
## Residual deviance: 723.45  on 759  degrees of freedom
## AIC: 741.45
## 
## Number of Fisher Scoring iterations: 5

Metode Stepwise: Forward, Backward, dan Kedua Arah

null_model <- glm(Outcome ~ 1, data = data, family = binomial)
step_model <- step(null_model, scope = formula(model_full), direction = "both", trace = FALSE)
summary(step_model)
## 
## Call:
## glm(formula = Outcome ~ Glucose + BMI + Pregnancies + DiabetesPedigreeFunction + 
##     BloodPressure + Age + Insulin, family = binomial, data = data)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -8.4051362  0.7167033 -11.727  < 2e-16 ***
## Glucose                   0.0351123  0.0036625   9.587  < 2e-16 ***
## BMI                       0.0900886  0.0144619   6.229 4.68e-10 ***
## Pregnancies               0.1231724  0.0320688   3.841 0.000123 ***
## DiabetesPedigreeFunction  0.9475954  0.2980063   3.180 0.001474 ** 
## BloodPressure            -0.0132136  0.0051537  -2.564 0.010350 *  
## Age                       0.0147888  0.0092897   1.592 0.111393    
## Insulin                  -0.0011570  0.0008142  -1.421 0.155275    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 993.48  on 767  degrees of freedom
## Residual deviance: 723.45  on 760  degrees of freedom
## AIC: 739.45
## 
## Number of Fisher Scoring iterations: 5

Interpretasi

Model hasil stepwise memilih variabel yang memberikan kontribusi signifikan terhadap Outcome berdasarkan kriteria AIC. Model akhir mencakup variabel: Glucose, BMI, Pregnancies, DiabetesPedigreeFunction, BloodPressure, Age, dan Insulin. Variabel SkinThickness dikeluarkan karena tidak menambah keakuratan model secara signifikan.

Evaluasi Model: ROC dan AUC

pred_prob <- predict(step_model, type = "response")
roc_obj <- pROC::roc(response = data$Outcome, predictor = pred_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(roc_obj, main = "Kurva ROC", col = "blue")

pROC::auc(roc_obj)
## Area under the curve: 0.8396

Interpretasi

Kurva ROC mengilustrasikan trade-off antara sensitivitas dan 1-spesifisitas untuk berbagai nilai ambang (threshold). Area Under Curve (AUC) digunakan sebagai ukuran diskriminasi model.

AUC = 0.8396, yang mengindikasikan bahwa model memiliki kinerja klasifikasi yang baik.

Pseudo R-Squared

DescTools::PseudoR2(step_model, which = "all")
##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.2718016       0.2556966       0.2964400       0.4084782       0.2601373 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.4612332       0.3278732       0.4405103       0.3259582     739.4534296 
##             BIC          logLik         logLik0              G2 
##     776.6037474    -361.7267148    -496.7419551     270.0304806
pscl::pR2(step_model)
## fitting null model for pseudo-r2
##          llh      llhNull           G2     McFadden         r2ML         r2CU 
## -361.7267148 -496.7419551  270.0304806    0.2718016    0.2964400    0.4084782
rcompanion::nagelkerke(step_model)
## $Models
##                                                                                                                               
## Model: "glm, Outcome ~ Glucose + BMI + Pregnancies + DiabetesPedigreeFunction + BloodPressure + Age + Insulin, binomial, data"
## Null:  "glm, Outcome ~ 1, binomial, data"                                                                                     
## 
## $Pseudo.R.squared.for.model.vs.null
##                              Pseudo.R.squared
## McFadden                             0.271802
## Cox and Snell (ML)                   0.296440
## Nagelkerke (Cragg and Uhler)         0.408478
## 
## $Likelihood.ratio.test
##  Df.diff LogLik.diff  Chisq    p.value
##       -7     -135.02 270.03 1.4999e-54
## 
## $Number.of.observations
##           
## Model: 768
## Null:  768
## 
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
## 
## $Warnings
## [1] "None"

Interpretasi

Beberapa nilai pseudo R-squared diperoleh: - McFadden: 0.27 → menunjukkan kecocokan model sedang hingga baik. - Cox & Snell: 0.296 - Nagelkerke: 0.408 → menunjukkan sekitar 41% variasi Outcome dijelaskan oleh model. - VeallZimmermann dan McKelveyZavoina juga menunjukkan nilai > 0.4, memperkuat bahwa model cukup representatif terhadap data.

Tabel Klasifikasi dan Evaluasi

pred_class <- ifelse(pred_prob >= 0.5, 1, 0)
conf_matrix <- caret::confusionMatrix(factor(pred_class), data$Outcome, positive = "1")
conf_matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 445 111
##          1  55 157
##                                          
##                Accuracy : 0.7839         
##                  95% CI : (0.753, 0.8125)
##     No Information Rate : 0.651          
##     P-Value [Acc > NIR] : 7.051e-16      
##                                          
##                   Kappa : 0.5001         
##                                          
##  Mcnemar's Test P-Value : 1.965e-05      
##                                          
##             Sensitivity : 0.5858         
##             Specificity : 0.8900         
##          Pos Pred Value : 0.7406         
##          Neg Pred Value : 0.8004         
##              Prevalence : 0.3490         
##          Detection Rate : 0.2044         
##    Detection Prevalence : 0.2760         
##       Balanced Accuracy : 0.7379         
##                                          
##        'Positive' Class : 1              
## 
conf_matrix$byClass[c("Sensitivity", "Specificity")]
## Sensitivity Specificity 
##   0.5858209   0.8900000

Interpretasi

  • Akurasi: 78.39%, berarti sekitar 78% observasi berhasil diprediksi dengan benar.
  • Sensitivity: 0.586 → model dapat mengenali 58.6% dari total kasus positif (penderita diabetes).
  • Specificity: 0.89 → model mampu mengenali 89% dari total kasus negatif.
  • Kappa: 0.5 menunjukkan kesesuaian sedang.
  • McNemar’s Test signifikan (p < 0.05) mengindikasikan ketidakseimbangan klasifikasi, terutama pada kesalahan klasifikasi positif.

Visualisasi Precision-Recall Curve (Opsional Jika Data Tidak Seimbang)

library(PRROC)
pr <- pr.curve(scores.class0 = pred_prob[data$Outcome == 1],
               scores.class1 = pred_prob[data$Outcome == 0],
               curve = TRUE)
plot(pr, main = "Precision-Recall Curve")

Kesimpulan Akhir

Model regresi logistik yang dibangun menunjukkan performa yang cukup baik dalam memprediksi status diabetes berdasarkan variabel-variabel medis. Dengan AUC ROC sebesar 0.8396 dan AUC PR sebesar 0.7288. Evaluasi klasifikasi menunjukkan akurasi sekitar 78%, dengan keseimbangan yang cukup antara sensitivitas dan spesifisitas. Pseudo R-squared juga mengindikasikan model cukup representatif. Model ini dapat digunakan sebagai alat bantu prediktif yang baik dalam konteks medis, khususnya untuk skrining risiko diabetes.

Referensi

  1. Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Wiley.

  2. Agresti, A. (2007). An Introduction to Categorical Data Analysis (2nd ed.). Wiley-Interscience.

  3. Dobson, A. J., & Barnett, A. (2018). An Introduction to Generalized Linear Models (4th ed.). CRC Press.

  4. McCullagh, P., & Nelder, J. A. (1989). Generalized Linear Models (2nd ed.). Chapman and Hall/CRC.

  5. Fox, J. (2015). Applied Regression Analysis and Generalized Linear Models (3rd ed.). Sage Publications.

  6. Cameron, A. C., & Trivedi, P. K. (2013). Regression Analysis of Count Data (2nd ed.). Cambridge University Press.

  7. Hilbe, J. M. (2011). Negative Binomial Regression (2nd ed.). Cambridge University Press.

  8. Menard, S. (2002). Applied Logistic Regression Analysis (2nd ed.). Sage Publications.

  9. Friendly, M., & Meyer, D. (2016). Discrete Data Analysis with R. CRC Press.

  10. Jaya, I. G. N. M. (2022). Analisis Data Kategori. Departemen Statistika FMIPA Universitas Padjadjaran.

  11. Subana, M., & Sudrajat, E. (2005). Dasar-Dasar Statistika Penelitian. Pustaka Setia.

  12. Sudjana, N. (2005). Metoda Statistika. Tarsito.

  13. Walpole, R. E., Myers, R. H., Myers, S. L., & Ye, K. (2012). Ilmu Peluang dan Statistik untuk Insinyur dan Ilmuwan (ed. Indonesia). Penerbit Erlangga.

  14. Hartono, J. (2005). Metode Penelitian Bisnis: Salah Kaprah dan Pengalaman-pengalaman. BPFE UGM.

  15. Hasan, M. I. (2002). Pokok-Pokok Materi Statistik 2 (Statistik Induktif). Bumi Aksara.

  16. O’Hara, R. B., & Kotze, D. J. (2010). Do not log-transform count data. Methods in Ecology and Evolution, 1(2), 118–122. https://doi.org/10.1111/j.2041-210X.2010.00021.x