Kata Pengantar

Assalamu’alaikum Warahmatullahi Wabarakatuh

Puji Syukur Kehadirat Allah SWT atas anugerah dan pertolongan-Nya sehingga penulis dapat menyelesaikan buu berjudul “Analisis Data Kategori”

Tak lupa juga penulis menyampaikan terima kasih sebesar-besarnya untuk Bapak Dr. I Gede Nyoman Mindra Jaya, M.Si selaku dosen telah memberi arahan selama penulis menyelesaikan buku ini.

Buku “Analisis Data Kategori” disusun sebagai bahan pembelajaran untuk memahami konsep dasar, metode, serta penerapan dari analisis data kategori yang sering ditermukan dalam berbagai bidang seperti ilmu sosial, kesehatan dan masih banyak lagi.

Penulis berharap buku ini bisa bermanfaat bagi pembaca dan menjadi sumber referensi yang berguna dalam memahami analisis data kategori. Namun, dalam pembuatan buku ini, peneliti menyadari bahwa buku ini tak lepas dari kekurangan. Oleh karena itu, penulis menyampaikan permohonan maaf serta terbuka untuk kritik dan saran demi perbaikan di masa mendatang.

Wassalamu’alaikum Warahmatullahi Wabarakatuh

Pendahuluan

Data kategorik merupakan jenis data yang terdiri dari nilai nilai yang dikelompokkan ke dalam kategori. Variabel kategorik memiliki skala pengukuran yang terdiri dari sekumpulan kategori. Data ini bisa berupa nominal ataupun ordinal.

Tujuan Analisis Data Kategori

Dapat digunakan untuk mengidentifikasi pola tertentu
Untuk mengetahui korelasi atau hubungan antar variabel
Sebagai alat bantu dalam pengambilan keputusan
Membantu pengembangkan model prediktif

Manfaat Analisis Data Kategori

Analisis data kategori memiliki manfaat dalam berbagai bidang, seperti bidang akademik, bidang ilmu sosial, bidang kesehatan, bidang dan masih banyak lagi. Contoh penerapannya dalam bidang akademik berupa survey kepuasan mahasiswa terhadap metode pengajaran dosen. Dengan adanya survey ini, dosen bisa mengevaluasi metode pembelajaran sehingga untuk tahun depan dosen bisa menggunakan metode pengajaran yang lebih baik lagi.

Metode dalam Analisis Data Kategori

Berbagai metode dapat digunakan dalam analisis data kategori, tergantung pada tujuan penelitian. Beberapa metode umum meliputi:

Tabel Kontingensi dan Uji Chi-Square

Digunakan untuk menguji hubungan natar dua variabel kategori
Contohnya, menguji apakah ada hubungan antara tingkat pendidikan dan status pekerjaan.

Regresi Logistik

Digunakan untuk memprediksi probabilitas suatu kejadian berdasarkan variabel kategori.
Contohnya, memprediksi apakah seorang pelanggan akan membeli produk berdasarkan kategori pref erensi mereka.

Analisis Korespondensi

Digunakan untuk memahami hubungan antara berbagai kategori dalam satu dataset.
Contohnya, analisis tentang preferensi makanan berdasarkan kelompok usia.

Decision Tree dan Random Forest

Metode machine learning yang sering digunakan untuk klasifikasi berbasis kategori.
Contohnya, mengklasifikasikan pelanggan berdasarkan tingkat loyalitas mereka.

Analisis data kategori merupakan bagian penting dari analisis statistik yang memiliki aplikasi luas dalam berbagai bidang. Dengan pendekatan yang tepat, analisis ini dapat membantu dalam memahami pola, hubungan, dan tren dalam data yang bersifat kategori. Selain itu, metode analisis data kategori seperti tabel kontingensi, uji chi-square, regresi logistik, dan machine learning terus berkembang seiring dengan kemajuan teknologi, memungkinkan pengambilan keputusan yang lebih akurat dan berbasis data. Dengan demikian, pemahaman yang baik tentang analisis data kategori sangat penting bagi peneliti, praktisi bisnis, pembuat kebijakan, dan profesional dari berbagai disiplin ilmu untuk menghasilkan wawasan yang lebih dalam dan keputusan yang lebih tepat.

Distribusi Probabilitas dalam Data Kategori

Distribusi Bernoulli

Distribusi Bernoulli biasanya digunakan untuk percobaan biner yaitu percobaan yang hanya memiliki dua kemungkinan hasil seperti Sukses (1) dengan probabilitas p dan Gagal (0) dengan probabilitas 1-p dengan fungsi probabilitas nya adalah :

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

dengan : X = Variabel acak biner (0 atau 1), p = Probabilitas sukses (X=1), 1-p = Probabilitas gagal

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, maka distribusi Binomial memiliki fungsi probabilitas:

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

dengan : X = Jumlah keberhasilan dalam n percobaan, n = Jumlah percobaan, k = Jumlah keberhasilan, p = Probabilitas keberhasilan dalam satu percobaan

Distribusi Multinomial

Distribusi Multinomial adalah generalisasi lebih lanjut dari distribusi Binomial, digunakan ketika setiap percobaan memiliki lebih dari dua kemungkinan hasil. Jika suatu eksperimen dilakukan n kali, dan setiap percobaan dapat menghasilkan salah satu dari k kategori dengan probabilitas p1,p2, … ,pk, maka distribusi probabilitasnya:

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

dengan : Xi = Frekuensi kemunculan kategori ke-i, n = Jumlah total percobaan, xi = Jumla kejadian kategori ke-i, pi = Probabilitas kategori ke-i,

Distribusi Poisson

Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam interval waktu atau ruang tertentu dengan rata-rata kejadian per unit waktu/ruang. Fungsi probabilitasnya:

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

dengan : X = Jumlah kejadian dalam interval tertentu, lambda = Rata-rata kejadian dalam interval tersebut, k = Jumlah kejadian yang diamati,

Desain Sampling dalam Analisis Data Kategori

Eksperimen

Eksperimen adalah suatu metode penelitian di amna peneliti secara aktif mengontrol ataupun memanipulasi satu atau lebih variabel independen untuk mengamati dampak nya terhadap variabel dependen. Teknik sampling yang digunakan meliputi :

Simple Random Sampling
Stratified Random Sampling
Cluster Sampling

Studi Cohort

Studi Cohort Prospective

Studi kohort prospective merupakan suatu metode penelitian dengan karakteristik tertentu yang diikuti dari waktu ke waktu untuk mengamati kejadian yang ingin dipelajari. Teknik sampling yang digunakan meliputi :

Census Sampling
Systematic Sampling
Matched Sampling

Studi Cohort Retrospective

Studi kohort retrospektif, data historis digunakan untuk mengelompokkan individu berdasarkan paparan dan kemudian menganalisis hasil yang terjadi. Teknik sampling yang digunakan meliputi :

Convenience Sampling
Quota Sampling
Case-Based Sampling

Studi Kasus Kontrol

Studi Kasus-Kontrol merupakan jenis studi observasional yang digunakan untuk menyelidiki hubungan anatar faktor risiko dan suatu outcome. Berbeda dengan studi kohort, studi ini dimulai dari outcome dan bekerja mundur untuk menentukan paparan sebelumnya. Teknik sampling yang diugnakan meliputi :

Purposive Sampling
Snowball Sampling
Incidence Density Sampling

Tabel Kontingensi 2 x 2

Tabel kontingensi 2 x 2 memiliki struktur sebagai berikut :

Data untuk Tabel Kontingensi 2 x 2

Data = Studi Membandingkan Efek Konsumsi Aspirin terhadap Kejadian Serangan Jantung

Sumber Data = Agresti, Alan. An Introduction to Categorical Data Analysis, Second Edition. London: Wiley, 2007. dan Preliminary report: Findings from the aspirin component of the ongoing Physicians’ Health Study

data2x2<- matrix(c(189, 10845, 104, 10933), nrow = 2, byrow = TRUE)
dimnames(data2x2) <- list("Kelompok" = c("Placebo", "Aspirin"), "Kejadian" = c("Serangan Jantung (Ya)", "Serangan Jantung (Tidak)"))
print(data2x2)

##          Kejadian
## Kelompok  Serangan Jantung (Ya) Serangan Jantung (Tidak)
##   Placebo                   189                    10845
##   Aspirin                   104                    10933

Distribusi Peluang dalam Tabel Kontingensi 2 x 2

Peluang Bersama

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

Peluang Marjinal

Peluang marginal adalah probabilitas kejadian suatu variabel tanpa mempertimbangkan variabel lainnya

Peluang Bersyarat

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

Perhitungan menggunakan R

n <- sum(data2x2)
# Peluang Bersama
P_joint <- data2x2 / n
# Peluang Marginal
P_marginal_rows <- rowSums(data2x2) / n
P_marginal_cols <- colSums(data2x2) / n
# Peluang Bersyarat
P_conditional <- data2x2 / rowSums(data2x2)
list(Peluang_Bersama = P_joint, Peluang_Marginal_Baris = P_marginal_rows,
     Peluang_Marginal_Kolom = P_marginal_cols,
     Peluang_Bersyarat = P_conditional)

## $Peluang_Bersama
##          Kejadian
## Kelompok  Serangan Jantung (Ya) Serangan Jantung (Tidak)
##   Placebo           0.008563273                0.4913688
##   Aspirin           0.004712066                0.4953559
## 
## $Peluang_Marginal_Baris
##  Placebo  Aspirin 
## 0.499932 0.500068 
## 
## $Peluang_Marginal_Kolom
##    Serangan Jantung (Ya) Serangan Jantung (Tidak) 
##               0.01327534               0.98672466 
## 
## $Peluang_Bersyarat
##          Kejadian
## Kelompok  Serangan Jantung (Ya) Serangan Jantung (Tidak)
##   Placebo            0.01712887                0.9828711
##   Aspirin            0.00942285                0.9905771

Interpretasi :

Peluang Bersama

Dapat disimpulkan bahwa 0.86% dari seluruh peserta berada di kelompok Placebo dan mengalami serangan jantung dan
Dapat disimpulkan bahwa 0.47% dari seluruh peserta berada di kelompok Aspirin dan mengalami serangan jantung.

Peluang Marginal Baris

Dapat disimpulkan bahwa dari peserta hampir terbagi rata antara kelompok kelompok yang ada yaitu Placebo dan Aspirin

Peluang Marginal Kolom

Dapat disimpulkan bahwa dari keseluruhan peserta hanya 1.3% yang mengalami serangan jantung

Peluang Bersyarat

Di antara peserta yang diberi Placebo, 1.71% mengalami serangan jantung
Di antara peserta yang diberi Aspirin, 0,94% mengalami serangan jantung
Dari sini dapat dilihat bahwa Aspirin menurunkan risiko terkena serangan jantung

Ukuran Asosiasi dalam Dat Kategori 2x2

Risk Difference

Risk Difference (RD) atau Perbedaan Risiko adalah ukuran absolut dalam epidemiologi yang menggambarkan perbedaan antara probabilitas kejadian suatu hasil dalam dua kelompok yang berbeda. RD dihitung sebagai selisih antara risiko kejadian dalam kelompok terpapar dan risiko kejadian dalam kelompok tidak terpapar.

\[ RD = \frac{n11}{n1.} - \frac{n21}{n2.} \]

Jika RD > 0, maka risiko kejadian lebih tinggi di Grup 1 dibandingkan Grup 2.
Jika RD < 0, maka risiko kejadian lebih rendah di Grup 1 dibandingkan Grup 2.
Jika RD = 0, maka tidak ada perbedaan risiko antara dua kelompok

Relative Risk

Relative Risk (RR) atau Risiko Relatif adalah ukuran yang digunakan dalam epidemiologi untuk membandingkan risiko kejadian suatu peristiwa (misalnya penyakit atau kondisi tertentu) antara dua kelompok, yaitu kelompok yang terpapar dan kelompok yang tidak terpapar.

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

Jika RR > 1, maka kejadian lebih sering terjadi di Grup 1 dibandingkan Grup 2.
Jika RR < 1, maka kejadian lebih jarang terjadi di Grup 1 dibandingkan Grup 2.
Jika RR = 1, maka tidak ada perbedaan risiko antara dua kelompok.

Odds Ratio

Odds Ratio (OR) atau Rasio Odds adalah ukuran yang digunakan dalam epidemiologi dan statistik untuk membandingkan odds (peluang) terjadinya suatu kejadian antara dua kelompok, yaitu kelompok yang terpapar dan kelompok yang tidak terpapar.

\[ OR = \frac{n11 \cdot n22}{n12 \cdot n21} \]

Jika OR > 1, maka peluang kejadian lebih besar di Grup 1 dibandingkan Grup 2.
Jika OR < 1, maka peluang kejadian lebih kecil di Grup 1 dibandingkan Grup 2.
Jika OR = 1, maka tidak ada perbedaan peluang kejadian antara dua kelompok.

Perhitungan Menggunakan R

n11 = 189
n12 = 10845
n21 = 104
n22 = 10933

RD <- (n11 / (n11 + n12))- (n21 / (n21 + n22))
RR <- (n11 / (n11 + n12)) / (n21 / (n21 + n22))
OR <- (n11 * n22) / (n12 * n21)
list(Risk_Difference = RD, 
     Relative_Risk = RR,
     Odds_Ratio = OR)

## $Risk_Difference
## [1] 0.007706024
## 
## $Relative_Risk
## [1] 1.817802
## 
## $Odds_Ratio
## [1] 1.832054

Interpretasi

Risk Difference

Perbedaan risiko serangan jantung antara kelompok Placebo dan Aspirin adalah 0.77% sehingga aspirin menurunkan risiko ablosut sebesar 0.77%

Relative Risk

Risiko serangan jantung pada kelompok Placebo adalah sekitar 1.82 kali lebih tinggi dibandingkan kelompok Aspirin

Odds Ratio

Peluang mengalami serangan jantung pada kelompok Placebo sekitar 1.83 kali lebih besar dibanding kelompok Aspirin.

Inferensi Tabel Kontingensi Dua Arah

Inferensi dalam tabel kontingensi dua arah, inferensi digunakan untuk emnganalisis hubungan anatar dua variabel kategorikal. Inferensi dalam tabel kontingensi dua arah ada dua yaitu Estimasi dan Pengujian Hipotesis

Estimasi

Estimasi berfungsi untuk memperkirakan parameter populasi berdasarkan data sampel yang ada. Estimasi ada dua meliputi :

Estimasi Titik

Estimasi titik digunakan untuk menentukan satu nilai spesifik sebagai perkiraan terbaik dari parameter populasi

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

Estimasi Interval

Dalam estimasi titik, kita perlu menentukan nilai spesifik untuk memperkirakan parameter populasi. Nah, dalam estimasi interval dalam memperkirakan paramater populasi melalui rentang nilai dengan tingkat kepercayaan tertentu

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

Perhitungan dengan R

# Data
n11 = 189
n12 = 10845
n21 = 104
n22 = 10933
n1. = n11 + n12
n2. = n21 + n22

# Estimasi titik proporsi
p_hat_placebo <- n11 / n1.
p_hat_aspirin <- n21/ n2.

# Z-score untuk 95% CI
z <- 1.96

# Standard Error
se_placebo <- sqrt(p_hat_placebo * (1 - p_hat_placebo) / n1.)
se_aspirin <- sqrt(p_hat_aspirin * (1 - p_hat_aspirin) / n2.)

# Confidence Interval (CI)
ci_placebo <- c(p_hat_placebo - z * se_placebo, p_hat_placebo + z * se_placebo)
ci_aspirin <- c(p_hat_aspirin - z * se_aspirin, p_hat_aspirin + z * se_aspirin)

# Output hasil
list(
  Placebo = list(
    Estimasi_Titik = p_hat_placebo,
    CI_95 = ci_placebo
  ),
  Aspirin = list(
    Estimasi_Titik = p_hat_aspirin,
    CI_95 = ci_aspirin
  )
)

## $Placebo
## $Placebo$Estimasi_Titik
## [1] 0.01712887
## 
## $Placebo$CI_95
## [1] 0.01470783 0.01954992
## 
## 
## $Aspirin
## $Aspirin$Estimasi_Titik
## [1] 0.00942285
## 
## $Aspirin$CI_95
## [1] 0.00762039 0.01122531

Pengujian Hipotesis

Uji Proporsi

Uji proporsi digunakan untuk membandingkan proporsi kejadian antara dua kelompok dalam tabel kontingensi, terutama untuk menentukan apakah terdapat perbedaan yang signifikan dalam proporsi kejadian antara dua kelompok yang berbeda.

Hipotesis :

H0 :Tidak terdapat perbedaan proporsi antara dua kelompok

H1 :Terdapat perbedaan proporsi antara dua kelompok

Proporsi Sampel

\[ \hat{p}_1 = \frac{n_{11}}{n_{1.}}, \quad \hat{p}_2 = \frac{n_{21}}{n_{2.}} \]

Proporsi Gabungan

\[ \hat{p} = \frac{n_{11} + n_{21}}{n_{1.} + n_{2.}} \]

Statistik Uji Z

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

Perhitungan menggunakan R

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

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(data2x2[1, 1], data2x2[2, 1]) out of c(sum(data2x2[1, ]), sum(data2x2[2, ]))
## X-squared = 24.429, df = 1, p-value = 7.71e-07
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.004597134 0.010814914
## sample estimates:
##     prop 1     prop 2 
## 0.01712887 0.00942285

Interpretasi :

Dikarekan p-value(7.71e-07) < 0,05, maka dapat disimpulkan bahwa terdapat perbedaan proporsi kejadian antara kelompok yang diberi placebo dan kelompok yang diberi aspirin

Uji Asosiasi

Hipotesis :

H0 :Tidak terdapat asosiasi antara dua variabel

H1 :terdapat asosiasi antara dua variabel

Risk Difference

\[ RD = \frac{n11}{n1.} - \frac{n21}{n2.} \]

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 RD

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

Relative Risk

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

Standard Error untuk RR

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

Statistik Uji Z RR

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

Odds ratio

\[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot 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 OR

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

Perhitungan Menggunakan R

# Data
n11 = 189
n12 = 10845
n21 = 104
n22 = 10933
n1. = n11 + n12
n2. = n21 + n22

# Risk Difference
p1<-(n11/n1.)
p2<-(n21/n2.)
rd <- p1- p2
se_rd <- sqrt((p1 * (1- p1) / n1.) + p2*((1- p2) / n2.))
z_rd <- rd / se_rd
 
# Relative Risk
rr <- (n11/n1.) / (n21/n2.)
se_ln_rr <- sqrt((1/n11)- (1/n1.) + (1/n21)- (1/n2.))
z_rr <- log(rr) / se_ln_rr

# Odds Ratio
or <- (n11 * n22) / (n12 * n21)
se_ln_or <- sqrt((1/n11) + (1/n12) + (1/n21) + (1/n22))
z_or <- log(or) / se_ln_or
 
# Hasil
list(RD = rd, SE_RD = se_rd, Z_RD = z_rd, RR = rr, SE_Ln_RR = se_ln_rr, Z_RR = z_rr, OR = or, SE_ln_OR = se_ln_or, Z_OR = z_or)

## $RD
## [1] 0.007706024
## 
## $SE_RD
## [1] 0.001539964
## 
## $Z_RD
## [1] 5.00403
## 
## $RR
## [1] 1.817802
## 
## $SE_Ln_RR
## [1] 0.1213473
## 
## $Z_RR
## [1] 4.92494
## 
## $OR
## [1] 1.832054
## 
## $SE_ln_OR
## [1] 0.1228416
## 
## $Z_OR
## [1] 4.928604

Interpretasi :

Dari ketiga uji hipotesis asosiasi, ketiganya memiliki hasil yang signifikan sehingga dapat disimpulkan bahwa terdapat asosiasi dari kedua variabel

Uji Independensi

Uji Chi-Square

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

Langkah-langkah :

Menghitung nilai ekspektasi

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

dengan : Ri = total baris ke-i Cj = total kolom k-j N = total sampel

Menghitung Chi-Square

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

dengan : O = nilai observasi E = nilai ekspektasi

Perhitungan menggunakan R

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

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data2x2
## X-squared = 24.429, df = 1, p-value = 7.71e-07

Interpretasi :

Dikarenakan p-value (7.71e-07) < 0,05, maka dapat disimpulkan bahwa terdpat hubungan antara variabel Kelompok dan Kejadian

Uji Likelihood Ratio

Uji Likelihood Ratio merupakan alternatif dari uji Chi-Square dalam menguji apakah terdapt hubungan antara kedua variabel

Langkah-Langkah : 1. Mencari nilai ekspektasi mirip dengan Chi-Square 2. Menghitung Statistik Uji Likelihodd Ratio

Rumus Uji Likelihood Ratio

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

dengan : n ij = frekuensi observasi myu ij = frekuensi ekspektasi

Perhitungan dengan R

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


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

# Nilai kritis chi-square untuk df = 1 dan alpha = 0.05
critical_value <- qchisq(0.95, df = 1)


# Hasil
list(G2 = G2, Critical_Value = critical_value)

## $G2
## [1] 25.37196
## 
## $Critical_Value
## [1] 3.841459

Interpretasi :

Dikarenakan p-value < 0,05, maka dapat disimpulkan bahwa terdapat hubungan antar variabel kelompok dan kejadian

Uji Exact Fisher

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

Perhitungan menggunakan R

fisher.test(data2x2)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  data2x2
## p-value = 5.033e-07
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.432396 2.353927
## sample estimates:
## odds ratio 
##   1.831993

Interpretasi :

Dikarenakan p-value < 0,05 maka dapat disimpulkan abhwa terdapat hubungan yang signifikan antara kedua variabel

Analisis Residual dalam Tabel Kontingensi

Residual dalam tabel kontingensi digunakan untuk mengidentifikasi sel mana yang menyumbang paling banyak terhadap hubungan antara variabel kategori.

Jika residual 0 berarti tidak terdapt perbedaan signifikan antara jumlah observasi yang terjadi di suatu sel dengan jumlah yang diprediksi
Jika residual positif besar menunjukkan bahwa terdapat hubungan positif yang kuat di antar kedua variabel kategori tersebut
Jika residual negatif besar dapat menunjukkan hubungan negatif ataupun bisa juga menunjukkan bahwa tidak adanya asosiasi antar kedua variabel

Jenis Residual

Pearson Residual

\[ e_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}}} \]

Standardized Residual

\[ r_{ij} = \frac{O_{ij} - E_{ij}}{\sqrt{E_{ij}(1 - p_{i+})(1 - p_{+j})}} \]

Perhitungan menggunakan R

# Misal kamu sudah punya observed (tabel kontingensi):
observed <- data2x2

# Menghitung expected frequency
expected <- outer(rowSums(observed), colSums(observed)) / sum(observed)

# Pearson Residual
pearson_residual <- (observed - expected) / sqrt(expected)

# Standardized Residual
row_sum <- rowSums(observed)
col_sum <- colSums(observed)
total_sum <- sum(observed)

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

# Menampilkan hasil
list(
  Pearson_Residual = pearson_residual,
  Standardized_Residual = standardized_residual
)

## $Pearson_Residual
##          Kejadian
## Kelompok  Serangan Jantung (Ya) Serangan Jantung (Tidak)
##   Placebo              3.513202               -0.4075003
##   Aspirin             -3.512724                0.4074449
## 
## $Standardized_Residual
##          Kejadian
## Kelompok  Serangan Jantung (Ya) Serangan Jantung (Tidak)
##   Placebo              5.001388                -5.001388
##   Aspirin             -5.001388                 5.001388

Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah menganalisis hubungan antara tiga variabel kategori secara simultan. Dalam banyak situasi, hubungan antara dua variabel (misalnya X dan Y) dapat dipengaruhi oleh variabel ketiga 𝑍, yang disebut sebagai variabel kontrol atau kovariat.

Tabel Parsial dan Marginal

Struktur Tabel Kontingensi Tiga Arah

Data untuk Tabel Kontingensi 3 arah

Data : Perbandingan lamanya pendidikan (tahun) dan jenis kelamin dalam menentukan sikap terhadap peran wanita

Sumber Data : Maryana. (2013). Model Log Linier yang Terbaik untuk Analisis Data Kualitatif pada Tabel Kontingensi Tiga Arah. Malikussaleh Industrial Engineering Journal, 2(2), 32–37.

Keterangan :
Pernyataan yang diberikan “Wanita sebaiknya menjaga atau memelihara rumah mereka danikut suami selamanya”

library(knitr)
library(Epi)

## Warning: package 'Epi' was built under R version 4.4.3

data3arah <- array(
  c(72, 47,  
    158, 85,  
    110, 196,
    173, 283,   
    44, 179,
    28, 187),  
  dim = c(2, 2, 3), 
  dimnames = list(
    Jenis_Kelamin = c("Pria", "Wanita"),
    Jawaban = c("Setuju", "Tidak"),
    Pendidikan = c("<=8", "9-12", ">=13")
  )
)
data3arah

## , , Pendidikan = <=8
## 
##              Jawaban
## Jenis_Kelamin Setuju Tidak
##        Pria       72   158
##        Wanita     47    85
## 
## , , Pendidikan = 9-12
## 
##              Jawaban
## Jenis_Kelamin Setuju Tidak
##        Pria      110   173
##        Wanita    196   283
## 
## , , Pendidikan = >=13
## 
##              Jawaban
## Jenis_Kelamin Setuju Tidak
##        Pria       44    28
##        Wanita    179   187

Freq Parsial Berdasarkan Lama Pendidikan

Pendidikan <=8

freq_parsial_1 = data3arah[, , "<=8"] 
freq_parsial_1

##              Jawaban
## Jenis_Kelamin Setuju Tidak
##        Pria       72   158
##        Wanita     47    85

Pendidikan 9-12

freq_parsial_2 = data3arah[, , "9-12"] 
freq_parsial_2

##              Jawaban
## Jenis_Kelamin Setuju Tidak
##        Pria      110   173
##        Wanita    196   283

Pendidikan >=13

freq_parsial_3 = data3arah[, , ">=13"] 
freq_parsial_3

##              Jawaban
## Jenis_Kelamin Setuju Tidak
##        Pria       44    28
##        Wanita    179   187

Distribusi Peluang

Peluang Bersama

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

Perhitungan menggunakan R

# Hitung total keseluruhan
total <- sum(data3arah)

# Hitung peluang bersama
peluang_bersama <- data3arah / total
peluang_bersama

## , , Pendidikan = <=8
## 
##              Jawaban
## Jenis_Kelamin     Setuju      Tidak
##        Pria   0.04609475 0.10115237
##        Wanita 0.03008963 0.05441741
## 
## , , Pendidikan = 9-12
## 
##              Jawaban
## Jenis_Kelamin     Setuju     Tidak
##        Pria   0.07042254 0.1107554
##        Wanita 0.12548015 0.1811780
## 
## , , Pendidikan = >=13
## 
##              Jawaban
## Jenis_Kelamin     Setuju      Tidak
##        Pria   0.02816901 0.01792574
##        Wanita 0.11459667 0.11971831

Peluang Marginal

marginal_X <- apply(peluang_bersama, 1, sum)
marginal_Z <- apply(peluang_bersama, 3, sum)

list(Marginal_X = marginal_X, Marginal_Z = marginal_Z)

## $Marginal_X
##      Pria    Wanita 
## 0.3745198 0.6254802 
## 
## $Marginal_Z
##       <=8      9-12      >=13 
## 0.2317542 0.4878361 0.2804097

Peluang Bersyarat

#(Pria, Pendidikan <=8)
P_XY <- sum(data3arah["Pria", , "<=8"]) / sum(data3arah)  # P(Jenis_Kelamin = Pria, Pendidikan = <=8)
P_ZXY <- data3arah["Pria", "Setuju", "<=8"] / sum(data3arah)  # P(Setuju, Pria, <=8)
P_Setuju_given_Pria_8 <- P_ZXY / P_XY

#(Wanita, Pendidikan <=8)
P_XY <- sum(data3arah["Wanita", , "<=8"]) / sum(data3arah)
P_ZXY <- data3arah["Wanita", "Setuju", "<=8"] / sum(data3arah)
P_Setuju_given_Wanita_8 <- P_ZXY / P_XY

#(Pria, Pendidikan 9-12)
P_XY <- sum(data3arah["Pria", , "9-12"]) / sum(data3arah)
P_ZXY <- data3arah["Pria", "Setuju", "9-12"] / sum(data3arah)
P_Setuju_given_Pria_912 <- P_ZXY / P_XY

#(Wanita, Pendidikan 9-12)
P_XY <- sum(data3arah["Wanita", , "9-12"]) / sum(data3arah)
P_ZXY <- data3arah["Wanita", "Setuju", "9-12"] / sum(data3arah)
P_Setuju_given_Wanita_912 <- P_ZXY / P_XY

#(Pria, Pendidikan >=13)
P_XY <- sum(data3arah["Pria", , ">=13"]) / sum(data3arah)
P_ZXY <- data3arah["Pria", "Setuju", ">=13"] / sum(data3arah)
P_Setuju_given_Pria_13 <- P_ZXY / P_XY

#(Wanita, Pendidikan >=13)
P_XY <- sum(data3arah["Wanita", , ">=13"]) / sum(data3arah)
P_ZXY <- data3arah["Wanita", "Setuju", ">=13"] / sum(data3arah)
P_Setuju_given_Wanita_13 <- P_ZXY / P_XY


list(P_Setuju_given_Pria_8,
     P_Setuju_given_Wanita_8,
     P_Setuju_given_Pria_912,
     P_Setuju_given_Wanita_912,
     P_Setuju_given_Pria_13,
     P_Setuju_given_Wanita_13)

## [[1]]
## [1] 0.3130435
## 
## [[2]]
## [1] 0.3560606
## 
## [[3]]
## [1] 0.3886926
## 
## [[4]]
## [1] 0.4091858
## 
## [[5]]
## [1] 0.6111111
## 
## [[6]]
## [1] 0.489071

Interpretasi :

Semakin tingginya pendidikan, semakin kecil kecenderungan untuk setuju tentang pernyataan

Tabel Peluang Bersyarat

data3arah <- array(
  c(72, 47,  
    158, 85,  
    110, 196,
    173, 283,   
    44, 179,
    28, 187),  
  dim = c(2, 2, 3), 
  dimnames = list(
    Jenis_Kelamin = c("Pria", "Wanita"),
    Jawaban = c("Setuju", "Tidak"),
    Pendidikan = c("<=8", "9-12", ">=13")))
 
ftable(data3arah)

##                       Pendidikan <=8 9-12 >=13
## Jenis_Kelamin Jawaban                         
## Pria          Setuju              72  110   44
##               Tidak              158  173   28
## Wanita        Setuju              47  196  179
##               Tidak               85  283  187

Ukuran Asosiasi

Risk Difference

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

Perhitungan menggunakan R

##Pendidikan <=8
p11 <- freq_parsial_1[1, 1] / sum(freq_parsial_1[1, ])
p21 <- freq_parsial_1[2, 1] / sum(freq_parsial_1[2, ])
RD1 <- p11 - p21


##Pendidikan 9-12
p12 <- freq_parsial_2[1, 1] / sum(freq_parsial_2[1, ])
p22<- freq_parsial_2[2, 1] / sum(freq_parsial_2[2, ])
RD2 <- p12 - p22

##Pendidikan >=13
p13 <- freq_parsial_3[1, 1] / sum(freq_parsial_3[1, ])
p23<- freq_parsial_3[2, 1] / sum(freq_parsial_3[2, ])
RD3 <- p13 - p23

list(RD_1 = RD1,RD_2 = RD2, RD_3 = RD3 )

## $RD_1
## [1] -0.04301713
## 
## $RD_2
## [1] -0.02049322
## 
## $RD_3
## [1] 0.1220401

Interpretasi :

Untuk lama pendidikan <=8 tahun kemungkinan pria untuk menjawab Setuju lebih rendah 4.3% dibanding wanita
Untuk lama pendidikan 9-12 tahun kemungkinan pria untuk menjawab Setuju lebih rendah 2% dibanding wanita
Untuk lama pendidikan <=8 tahun kemungkinan pria untuk menjawab Setuju lebih tinggi 12.2% dibanding wanita

Relative Risk

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

Perhitungan menggunakan R

##Pendidikan <=8 tahun
RR1 <- p11 / p21

##Pendidikan 9-12 tahun
RR2 <- p12 / p22

##Pendidikan >=13 tahun
RR3 <- p13 / p23


list(RR_1 = RR1, RR_2 = RR2, RR_3 = RR3)

## $RR_1
## [1] 0.8791859
## 
## $RR_2
## [1] 0.9499171
## 
## $RR_3
## [1] 1.249534

Interpretasi :

Untuk lama pendidikan <=8 tahun, pria memeiliki 12.1% lebih rendah kemungkinan untuk menjawab Setuju dibanding wanita
Untuk lama pendidikan 9-12 tahun, pria memeiliki 5% lebih rendah kemungkinan untuk menjawab Setuju dibanding wanita
Untuk lama pendidikan >=13 tahun, pria memeiliki 25% lebih tinggi kemungkinan untuk menjawab Setuju dibanding wanita

Odds Ratio

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

Perhitungan menggunakan R

##Pendidikan <=8 tahun
odds11 <- freq_parsial_1[1, 1] * freq_parsial_1[2, 2]
odds21 <- freq_parsial_1[1, 2] * freq_parsial_1[2, 1]
OR1 <- odds11 / odds21

##Pendidikan 9-12 tahun
odds12 <- freq_parsial_2[1, 1] * freq_parsial_1[2, 2]
odds22 <- freq_parsial_2[1, 2] * freq_parsial_1[2, 1]
OR2 <- odds12 / odds22

##Pendidikan >=13 tahun
odds13 <- freq_parsial_3[1, 1] * freq_parsial_3[2, 2]
odds23 <- freq_parsial_3[1, 2] * freq_parsial_3[2, 1]
OR3 <- odds13 / odds23

list(OR_1 = OR1, OR_2 = OR2, OR_3 = OR3)

## $OR_1
## [1] 0.8241314
## 
## $OR_2
## [1] 1.14992
## 
## $OR_3
## [1] 1.64166

Interpretasi :

- Untuk lama pendidikan <=8, Odds (kemungkinan relatif) pria untuk Setuju lebih kecil 17.6% dibanding wanita

- Untuk lama pendidikan 9-12, Odds (kemungkinan relatif) pria untuk Setuju lebih besar 15% dibanding wanita

- Untuk lama pendidikan >=13, Odds (kemungkinan relatif) pria untuk Setuju lebih besar 64.2% dibanding wanita

Inferensi

Uji CMH

Uji Cochran-Mantel-Haenszel (CMH) digunakan untuk menguji hubungan antara dua variabel kategori dengan mempertimbangkan efek dari variabel perancu (confounder). Uji ini berfungsi dalam :

Menguji hubungan yang dikendalikan oleh faktor perancu dalam tabel kontingensi berlapis (stratified contingency tables).
Menguji hipotesis independensi antara dua variabel kategori dengan mempertimbangkan efek dari variabel ketiga.
Mengatasi bias akibat confounding, misalnya dalam studi epidemiologi atau eksperimen sosial.

Hipotesis :

H0 : 𝜃𝑥𝑦(𝑘) = 1 untuk setiap 𝑘 = 1,2,…,K

H1 : 𝜃𝑥𝑦(𝑘) ≠ 1 untuk paling sedikit satu k

Rumus Statistik Uji CMH :

\[ CMH = \frac{\left[ \sum_k (n_{11k} - \mu_{11k}) \right]^2}{\sum_k \text{var}(n_{11k})} \]

Perhitungan menggunakan R

cmh_base <- mantelhaen.test(data3arah, correct = FALSE)
cmh_base

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  data3arah
## Mantel-Haenszel X-squared = 0.00010377, df = 1, p-value = 0.9919
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.800176 1.252638
## sample estimates:
## common odds ratio 
##          1.001165

Interpretasi :

Dikarenakan p-value > 0.05, maka dapat disimpulkan bahwa tidak terdapat hubungan yang signifikan antara jenis kelamin dan jawaban setuju/tidak nya. Perbedaan jawaban antara pria dan wanita dapat dijelaskan oleh faktor pendidikan

Uji Homogenitas Breslow-Day

Hipotesis :

H0 : Odds ratio sama di seluruh strata

H1 : Setidaknya ada satu odds ratio yang berbeda

Rumus Statistik Uji Breslow-Day

\[ X^2_{\text{HBD}} = \sum_{j=1}^{K} \frac{(a_j - \tilde{a}_j)^2}{\widehat{\text{Var}}(a_j \mid \widehat{OR}_{MH})} \]

Perhitungan menggunakan R :

# Buat list 2x2 tables dari array 3 dimensi
tabel_strata <- lapply(1:dim(data3arah)[3], function(i) data3arah[,,i])
names(tabel_strata) <- dimnames(data3arah)$Pendidikan

# Gabung semua 2x2 table jadi 1 array 3 dimensi
array_strata <- array(unlist(tabel_strata), dim = c(2, 2, 3))

# Cek dimensi
dimnames(array_strata) <- list(
  Jenis_Kelamin = c("Pria", "Wanita"),
  Jawaban = c("Setuju", "Tidak"),
  Pendidikan = names(tabel_strata)
)

# Lakukan Uji Breslow-Day
library(DescTools)

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

breslow_test <- BreslowDayTest(array_strata)

# Tampilkan hasil
print(breslow_test)

## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  array_strata
## X-squared = 4.5997, df = 2, p-value = 0.1003

Interpretasi :

Dikarenakan p-value > 0.05, maka tidak terdapat perbedaan singifikan pada odds ratio antar kelompok pendidikan.

Generalized Linear Model

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

Distribusi dari exponential family
Fungsi link
Fungsi linear prediktor

Exponential Family

Berikut distribusi exponential family

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

Regresi Logistik

Persamaan regresi logistik menyerupai regresi linear, di mana nilai input dikombinasikan secara linear dengan koefisien (bobot) untuk menghasilkan prediksi. Namun, regresi logistik membatasi hasil prediksi menjadi nilai biner, yaitu 0 atau 1, dengan menggunakan fungsi aktivasi sigmoid.

Fungsi Sigmoid :

\[ f(x) = \frac{1}{1 + e^{-x}} \]

Output yang diprediksi oleh model terletak dalam rentang antara 0 hingga 1 dan mengikuti bentuk kurva S (S-shaped).

Data

Data : Heart Disease

Sumber : UCI Machine Learning Repository

Tujuan : Untuk memprediksi apakah seseorang berisiko terkena penyakit jantung berdasarkan variabel yang ada

Variabel :

age = Umur
chol = kadar kolestrol
thalach = detak jantung maksimal saat treadmill

Perhitungan dengan R

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.4.3

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

## Warning: package 'tibble' was built under R version 4.4.3

## Warning: package 'tidyr' was built under R version 4.4.3

## Warning: package 'readr' was built under R version 4.4.3

## Warning: package 'purrr' was built under R version 4.4.3

## Warning: package 'dplyr' was built under R version 4.4.3

## Warning: package 'forcats' was built under R version 4.4.3

## Warning: package 'lubridate' was built under R version 4.4.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(caret)

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

## 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(readr)

# Baca file CSV dari direktori lokal (pastikan file ada di folder kerja R)
datareglog <- read.csv("dataregresilog.csv", header = TRUE, sep = ";", dec = ",")

datareglog$target <- sample(c(0,1), size = nrow(datareglog), replace = TRUE)

datareglog <- datareglog %>%
  mutate(target = factor(target, levels = c(0,1), labels = c("Tidak", "Ya")))

model <- glm(target ~ age + chol + thalach, data = datareglog, family = binomial)
summary(model)

## 
## Call:
## glm(formula = target ~ age + chol + thalach, family = binomial, 
##     data = datareglog)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)  1.002740   1.362585   0.736    0.462
## age         -0.018415   0.014327  -1.285    0.199
## chol         0.001596   0.002295   0.696    0.487
## thalach     -0.002323   0.005516  -0.421    0.674
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 419.89  on 302  degrees of freedom
## Residual deviance: 418.04  on 299  degrees of freedom
## AIC: 426.04
## 
## Number of Fisher Scoring iterations: 3

Interpretasi Odds Ratio

exp(coef(model))

## (Intercept)         age        chol     thalach 
##   2.7257402   0.9817535   1.0015974   0.9976798

Interpretasi :

Untuk variabel umur, setiap kenaikan 1 tahun usia, odds menderita penyakit akan turun sebesar 1,7%
Untuk variabel chol, setiap kenaikan 1 unit kolestrol, odds menderita penyakit kolestrol akan meningkat 0,08%
Untuk variabel thalach, setiap kenaikan 1 detak jantung maksimal saat treadmill, odds untuk terkena penyakit jantung turun sebesar 1,04%

Prediksi

# Prediksi probabilitas
datareglog$prob <- predict(model, type = "response")

# Konversi ke kelas berdasarkan threshold 0.5
datareglog$pred <- ifelse(datareglog$prob > 0.5, "Ya", "Tidak")
datareglog$pred <- factor(datareglog$pred, levels = c("Tidak", "Ya"))

# Evaluasi dengan confusion matrix
library(caret)
confusionMatrix(datareglog$pred, datareglog$target)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Tidak  Ya
##      Tidak    62  53
##      Ya       86 102
##                                           
##                Accuracy : 0.5413          
##                  95% CI : (0.4833, 0.5984)
##     No Information Rate : 0.5116          
##     P-Value [Acc > NIR] : 0.164321        
##                                           
##                   Kappa : 0.0774          
##                                           
##  Mcnemar's Test P-Value : 0.006644        
##                                           
##             Sensitivity : 0.4189          
##             Specificity : 0.6581          
##          Pos Pred Value : 0.5391          
##          Neg Pred Value : 0.5426          
##              Prevalence : 0.4884          
##          Detection Rate : 0.2046          
##    Detection Prevalence : 0.3795          
##       Balanced Accuracy : 0.5385          
##                                           
##        'Positive' Class : Tidak           
##

Regresi Poisson

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

Data

Sumber : Table 3.5 Reaven and Miller (1979) (Wiley Series in Probability and Statistics) Method of Multivariate Analysis by Alvin C

relative_weight (prediktor)
fasting_plasma_glucose (prediktor)
insulin_resistance (respon)

Perhitungan Menggunakan R

# Baca data dan ubah kolom menjadi terpisah
dataregpoi <- read.csv("datapasien.csv", sep = ";", header = TRUE, dec=",")
head(dataregpoi)

##   relative_weight fasting_plasma_glucose insulin_resistance
## 1            0.81                     80                 55
## 2            0.95                     97                 76
## 3            0.94                    105                105
## 4            1.04                     90                108
## 5            1.00                     90                143
## 6            0.76                     86                165

### Summary Data
summary(dataregpoi)

##  relative_weight  fasting_plasma_glucose insulin_resistance
##  Min.   :0.7200   Min.   : 66.00         Min.   : 32.00    
##  1st Qu.:0.8375   1st Qu.: 86.00         1st Qu.: 64.50    
##  Median :0.9450   Median : 90.00         Median : 92.00    
##  Mean   :0.9305   Mean   : 91.05         Mean   : 98.72    
##  3rd Qu.:1.0000   3rd Qu.: 98.00         3rd Qu.:128.00    
##  Max.   :1.2000   Max.   :106.00         Max.   :235.00

# Model Poisson
model_poisson <- glm(insulin_resistance ~ relative_weight + fasting_plasma_glucose,
                     data = dataregpoi, family = poisson())

summary(model_poisson)

## 
## Call:
## glm(formula = insulin_resistance ~ relative_weight + fasting_plasma_glucose, 
##     family = poisson(), data = dataregpoi)
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             4.845005   0.185254  26.153   <2e-16 ***
## relative_weight         1.293960   0.128260  10.089   <2e-16 ***
## fasting_plasma_glucose -0.016187   0.001882  -8.602   <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: 778.44  on 39  degrees of freedom
## Residual deviance: 639.23  on 37  degrees of freedom
## AIC: 898.44
## 
## Number of Fisher Scoring iterations: 4

Interpretasi :

Seluruh variabel prediktor sangat signifikan dalam model ini (p-value < 0,05)
Variabel berat relatif meningkatkan resiko resistensi insulin
Variabel gula darah setelah puasa menurunkan resiko resistensi insulin

 exp(coef(model_poisson))

##            (Intercept)        relative_weight fasting_plasma_glucose 
##            127.1039663              3.6472022              0.9839437

Interpretasi :

Untuk variabel berat relatif, setiap kenaikan satu unit laju resistensi insulin naik 3,65 kali
Untuk variabel gula darah setelah puasa, setiap kenaikan 1 unit gula darah setelah puasa laju resistensi insulin menurut sebesar 1,6%

Inferensi Generalized Linear Model

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

Ekspektasi dan Varians GLM

Ekspektasi

Jika diturunkan berdasarkan fungsi momen :

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

Untuk exponential family :

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

Maka turunan pertama :

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

Untuk ekspektasi turunan pertama :

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

Maka :

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

Untuk turunan kedua :

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

Sehingga :

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

Metode Penaksiran Parameter

Maximum Likelihood Estimation (MLE) Prinsip dasar : memaksimumkan fungsi likelihood/log-likelihood. Langkah:

Turunan pertama = 0
Turunan kedua < 0

Namun karena bentuk GLM tidak eksplisit, digunakan metode numerik.

Metode Optimisasi

Newton-Raphson

Menggunakan score vector (gradien)
Menggunakan Hessian matrix

Implementasi Newton-Raphton

Statistik score ke-j :

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

Turunan kedua :

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

Taylor Expansion

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

Estmasi Parameter :

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

Diagnostik Model GLM

Diagnostik digunakan untuk mengevaluasi apakah model sudah tepat dengan menggunakan

Uji formal
Grafik antar nilai prediksi dibandingkan dengan nilai aktual

Detail Metode Estimasi dan Inferensi Regresi Logistik

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

Fungsi Model Logistik

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

Log-likelihood untuk n observasi

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

Langkah-Langkah Newton-Raphson

Turunan Pertama (Score Function)\[ \mathbf{U}(\beta) = \frac{\partial \ell(\beta)}{\partial \beta} = \mathbf{X}^\top (\mathbf{y} - \boldsymbol{\pi}) \]
Turunan Kedua (Hessian Matrix)

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

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

Estimasi MLE dengan Newton-Raphson

# Estimasi MLE dengan Newton-Raphson (Manual) 
set.seed(42)

# Data simulasi baru
n <- 150
x1 <- rnorm(n, mean = 3, sd = 2)
x2 <- runif(n, min = -1, max = 1)
X <- cbind(1, x1, x2)

beta_true <- c(0.5, -2, 1.5)
eta <- X %*% beta_true
p <- 1 / (1 + exp(-eta))
y <- rbinom(n, 1, p)

# Newton-Raphson Iterasi Manual
beta <- rep(0, 3)  # inisialisasi
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
  
  # Cek konvergensi
  if (sum(abs(beta_new - beta)) < tol) {
    cat("Konvergen pada iterasi ke-", i, "\n")
    break
  }
  beta <- beta_new
}

## Konvergen pada iterasi ke- 9

# Hasil estimasi
print(beta)

##           [,1]
##     0.03495636
## x1 -1.68847181
## x2  2.69421377

Uji Wald

Tujuan Uji Wald :

Untuk menguji signifikansi parameter $\hat{\beta}_j$ dalam model regresi logistik :

H0 : $\hat{\beta}_j$ = 0

Dari teori estimasi MLE, estimator $\hat{\beta}_j$ mendekati distribusi normal:

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

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

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

Dengan statistik Wald:

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

set.seed(123)
 n <- 100
 x <- rnorm(n)
 log_odds <--0.5 + 1.2 * x
 p <- 1 / (1 + exp(-log_odds))
 y <- rbinom(n, 1, p)
 data <- data.frame(x, y)

 model <- glm(y ~ x, data = data, family = binomial)
 summary(model)

## 
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3097     0.2296  -1.349    0.177    
## x             1.2663     0.3080   4.111 3.94e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 137.99  on 99  degrees of freedom
## Residual deviance: 114.76  on 98  degrees of freedom
## AIC: 118.76
## 
## Number of Fisher Scoring iterations: 4

Uji Likelihood Ratio

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

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

## Analysis of Deviance Table
## 
## Model 1: y ~ 1
## Model 2: y ~ x
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        99     137.99                          
## 2        98     114.76  1   23.229 1.438e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Detail Metode Estimasi dan Inferensi Regresi Poisson

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

Fungsi distribusi Poisson:

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

Model regresi Poisson:

\[ \log(\lambda_i) = \mathbf{x}_i^\top \boldsymbol{\beta} \]

Estimasi Parameter (MLE)

Log-likelihood fungsi:

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

Dengan:

\[ \lambda_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta}) \]

Estimasi dilakukan dengan metode iterasi (IRLS).

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

Tahap 1: Definisikan Model Regresi Poisson

\[ \log(\lambda_i) = \mathbf{x}_i^\top \boldsymbol{\beta} \quad \text{sehingga} \quad \lambda_i = \exp(\mathbf{x}_i^\top \boldsymbol{\beta}) \]

Tahap 2: Mencari Log-likelihood yang Dimaksimalkan

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

Tahap 3: Formulasi Iteratif

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

Dengan:

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

dan

\[ \eta_i = \log(\lambda_i) = \mathbf{x}_i^\top \boldsymbol{\beta} \]

Simulasi

set.seed(123)
 n <- 100
 x <- rnorm(n)
 X <- cbind(1, x) # Tambah intercept
 beta_true <- c(0.5, 0.8)
 eta <- X %*% beta_true
 lambda <- exp(eta)
 y <- rpois(n, lambda)

IRLS Manual

# Inisialisasi
 beta <- c(0, 0)
 tol <- 1e-6
 max_iter <- 100
 for (i in 1:max_iter) {
 eta <- X %*% beta
 lambda <- exp(eta)
 W <- diag(as.numeric(lambda))
 z <- eta + (y- lambda) / lambda
 beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
 if (sum(abs(beta_new- beta)) < tol) {
 cat("Konvergen pada iterasi ke-", i, "\n")
 break
 }
 beta <- beta_new
 }

## Konvergen pada iterasi ke- 8

# hasil estimasi
beta

##        [,1]
##   0.4494951
## x 0.8600048

#Perbandingan dengan glm()
model_glm <- glm(y ~ x, family = poisson)
summary(model_glm)

## 
## Call:
## glm(formula = y ~ x, family = poisson)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.44950    0.08872   5.066 4.05e-07 ***
## x            0.86000    0.07463  11.523  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.05  on 99  degrees of freedom
## Residual deviance: 106.78  on 98  degrees of freedom
## AIC: 325.76
## 
## Number of Fisher Scoring iterations: 5

Pengujian Hipotesis Uji Wald

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

## Z: 4.110965 
## Chi-Square: 16.90003 
## p-value: 3.940095e-05

Evaluasi Model (AIC & BIC)

AIC(model)

## [1] 118.7598

BIC(model)

## [1] 123.9701

Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio

Regresi logistik merupakan salah satu metode analisis statistik yang digunakan untuk memodelkan hubungan antara satu variabel respons biner (dua kategori) dengan satu atau lebih variabel prediktor. Dalam penerapannya, prediktor yang digunakan dapat memiliki skala pengukuran berbeda, yaitu:

Nominal: Variabel yang tidak memiliki urutan logis antar kategorinya, hanya sebagai pembeda. Contoh: jenis kelamin (laki-laki, perempuan). Dalam analisis regresi logistik, data nominal diubah menjadi variabel dummy (misal: Female = 1, Male = 0).
Ordinal: Variabel yang memiliki urutan logis antar kategori, tetapi jarak antar kategori belum tentu sama. Contoh: tingkat pendidikan (SMA < Sarjana < Master < Doktor). Dalam analisis, variabel ordinal bisa diperlakukan:
- Sebagai nominal dengan membuat dummy untuk setiap kategori.– Atau diperlakukan seperti rasio dengan memberi nilai tertentu (misalnya SMA = 1, Sarjana = 2, Master = 3, Doktor = 4) untuk mencerminkan urutan dan jarak antar tingkatannya.
Rasio: Variabel numerik kontinu yang memiliki nol absolut dan rasio bermakna. Contoh: pendapatan (dalam juta rupiah per bulan). Data ini dapat langsung digunakan dalam model tanpa transformasi khusus.

Simulasi Data

# Set seed untuk konsistensi
set.seed(164)

# Jumlah sampel
n <- 500

# Variabel-variabel
income <- rnorm(n, mean = 70000, sd = 20000)          
age <- round(rnorm(n, mean = 35, sd = 10))
age <- pmax(age, 20)
environmental_awareness <- sample(c("Low", "Medium", "High"), n, replace = TRUE, prob = c(0.3, 0.4, 0.3))  

# Model logit: peluang memilih mobil listrik
logit_p <- -6 + 
  0.00007 * income -     
  0.02 * age +           
  1.5 * (environmental_awareness == "High") + 
  0.7 * (environmental_awareness == "Medium")  

# Ubah ke probabilitas
p <- 1 / (1 + exp(-logit_p))

# Variabel respon: memilih mobil listrik (1) atau tidak (0)
choose_ev <- rbinom(n, 1, p)

# Buat data frame
library(tibble)
dataev <- tibble(choose_ev, income, age, environmental_awareness)

# Lihat beberapa data
head(dataev)

## # A tibble: 6 × 4
##   choose_ev  income   age environmental_awareness
##       <int>   <dbl> <dbl> <chr>                  
## 1         1  70601.    34 High                   
## 2         1  94478.    30 Medium                 
## 3         1 104315.    28 Medium                 
## 4         0  52001.    20 High                   
## 5         0  81101.    26 Low                    
## 6         0  78464.    35 Medium

Latar Belakang Data

Seiring dengan perubahan iklim di dunia yang mengkhawatirkan, mobil listrik menjadi alternatif kendaraan yang menarik. Namun, keputusan individu dalam memiliki kendaraan listrik dipengaruhi oleh beberapa faktor, diantaranya penghasilan, umur, dan juga kesadaran akan lingkungan alam.

Melalui analisis, ini diharapkan memahami konsumen yang berpotensi membeli mobil lsitrik sehingga dapat membantu produsen dalam membuat kebijakan yang sesuai terkait mobil listrik ini.

Data yang disimulasikan berisi 500 observasi dengan variabel penghasilan, umur, dan tingkat kesadaran akan lingkungan alam.

Eksplorasi Data

library(dplyr)

dataev %>%
  group_by(choose_ev) %>%
  summarise(
    Jumlah = n(),
    Rata2_Income = mean(income)
  )

## # A tibble: 2 × 3
##   choose_ev Jumlah Rata2_Income
##       <int>  <int>        <dbl>
## 1         0    340       64292.
## 2         1    160       83756.

Perlakuan Variabel Ordinal

Treat sebagai Ordinal (Dummy)

dataev_nominal <- dataev %>%
mutate(
environmental_awareness = factor(environmental_awareness, levels = c("Low", "Medium", "High"))
)
model_nominal <- glm(choose_ev ~ income + age + environmental_awareness, data = dataev_nominal, family = binomial)
summary(model_nominal)

## 
## Call:
## glm(formula = choose_ev ~ income + age + environmental_awareness, 
##     family = binomial, data = dataev_nominal)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -7.287e+00  8.007e-01  -9.101  < 2e-16 ***
## income                         7.591e-05  8.290e-06   9.157  < 2e-16 ***
## age                           -5.076e-03  1.245e-02  -0.408   0.6836    
## environmental_awarenessMedium  9.393e-01  3.131e-01   3.000   0.0027 ** 
## environmental_awarenessHigh    1.931e+00  3.214e-01   6.009 1.87e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 626.87  on 499  degrees of freedom
## Residual deviance: 457.44  on 495  degrees of freedom
## AIC: 467.44
## 
## Number of Fisher Scoring iterations: 5

Odds Ratio

# Odds Ratio
odds_ratio <- exp(coef(model_nominal))
odds_ratio

##                   (Intercept)                        income 
##                  0.0006842836                  1.0000759159 
##                           age environmental_awarenessMedium 
##                  0.9949367698                  2.5582130595 
##   environmental_awarenessHigh 
##                  6.8971147720

Interpretasi

Interpretasi Koefisien

Intercept (-7,287)

Ini adalah log-odds dasar bagi individu dengan kesadaran lingkungan rendah, income = 0, dan age = 0, kondisi yang tidak realistis, tapi tetap menjadi baseline matematis.
Odds Ratio = 0,0006842836

income (0,00007591)

Signifikan (p < 2e-16) , artinya income berpengaruh signifikan terhadap keputusan memiliki kendaraan listrik atau tidak.
Odds Ratio =1,00008, artinya dalam setiap kenaikan $1 itu meningkatkan peluang membeli mobil listrik sebesar 0,00008%. Semakin besar penghasilan yang didapat akan meningkatkan peluang nya lebih besar juga

age (-0,005076)

TIdak signifikan (p = 0,6836), artinya age tidak berpengaruh signifikan terhadap keputusan memiliki kendaraan listrik atau tidak.
Odds Ratio = 0,995, artinya setiap peningkatan umur 1 tahun menurunkan peluang membeli mobil listrik sebesar 0,005%

environmental_awareness medium (0,9393)

Signifikan (p = 0,0027), artinya environmental_awareness medium berpengaruh signifikan terhadap keputusan memiliki kendaraan listrik atau tidak.
Odds Ratio = 2,56, peluang memilih EV yang memiliki kepedulian terhadap alam yang sedang 2,56 kali lebih tinggi dibanding yang memiliki kepedulian terhadap alam yang rendah

environmental_awareness high (1,931)

Signifikan (p = 1,87e-09), artinya environmental_awareness high berpengaruh signifikan terhadap keputusan memiliki kendaraan listrik atau tidak.
Odds Ratio = 6,9, peluang memilih EV yang memiliki kepedulian terhadap alam yang tinggi 6,9 kali lebih tinggi dibanding yang memiliki kepedulian terhadap alam yang rendah

Hasil Goodness of FIt

Null Deviance (626,87) merupakan deviansi dari model tanpa prediktor

Residual Deviance (457,44) merupakan deviansi dari model dengan prediktor. Dengan penurunan dari null deviance ke residual deviance menunjukkan model membawa informasi yang cukup baik

AIC (467,44) digunakan untuk membandingkan model.

Kesimpulan

Semakin tinggi kependulian akan lingkungan alam, semakin tinggi kemungkinan memilih kendaraan listrik
Income yang lebih tinggi akan meningkatkan kemungkin memilih kendaraan listrik tetapi tidak begitu signifikan
Umur tidak memiliki pengaruh yang singnifikan

Treat Variabel sebagai Rasio (Numeric Rank)

dataev_numeric <- dataev %>%
  mutate(
    environmental_awareness_numeric = case_when(
      environmental_awareness == "Low" ~ 1,
      environmental_awareness == "Medium" ~ 2,
      environmental_awareness == "High" ~ 3
    )
  )

# Model regresi logistik
model_numeric_ev <- glm(choose_ev ~ income + age + environmental_awareness_numeric, 
                        data = dataev_numeric, 
                        family = binomial)

# Summary model
summary(model_numeric_ev)

## 
## Call:
## glm(formula = choose_ev ~ income + age + environmental_awareness_numeric, 
##     family = binomial, data = dataev_numeric)
## 
## Coefficients:
##                                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     -8.268e+00  8.695e-01  -9.509  < 2e-16 ***
## income                           7.587e-05  8.279e-06   9.164  < 2e-16 ***
## age                             -5.043e-03  1.245e-02  -0.405    0.685    
## environmental_awareness_numeric  9.688e-01  1.582e-01   6.122 9.24e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 626.87  on 499  degrees of freedom
## Residual deviance: 457.45  on 496  degrees of freedom
## AIC: 465.45
## 
## Number of Fisher Scoring iterations: 5

Odds Ratio

odds_ratio <- exp(coef(model_numeric_ev))
odds_ratio

##                     (Intercept)                          income 
##                    0.0002565574                    1.0000758734 
##                             age environmental_awareness_numeric 
##                    0.9949698302                    2.6347043899

Interpretasi

Interpretasi Koefisien

Intercept (-8,268)

Signifikan (p < 2e-16), artinya Intercept berpengaruh signifikan terhadap keputusan memiliki kendaraan listrik atau tidak.
Odds Ratio = 0,00026

income (0,00007587)

Signifikan (p < 2e-16), artinya income berpengaruh signifikan terhadap keputusan memiliki kendaraan listrik atau tidak.
Odds Ratio = 1,00008, artinya dalam setiap kenaikan $1 itu meningkatkan peluang membeli mobil listrik sebesar 0,00008%. Semakin besar penghasilan yang didapat akan meningkatkan peluang nya lebih besar juga

age (-0,005043)

Tidak Signifikan (p = 0,685), artinya age tidak berpengaruh signifikan terhadap keputusan memiliki kendaraan listrik atau tidak.
Odds Ratio = 0,995, artinya setiap peningkatan umur 1 tahun menurunkan peluang membeli mobil listrik sebesar 0,005%

environmental_awareness_numeric (0,9688)

Signifikan (p = 9,24e-10), artinya environmental_awareness_numeric berpengaruh signifikan terhadap keputusan memiliki kendaraan listrik atau tidak.
Odds Ratio = 2,635, setiap naik satu level misalkan dari Low ke Medium, ini akan membuat peluang untuk memilih EV lebih besar 2,635 kali

Interpretasi Goodness of Fit

Null Deviance (626,87)

Residual Deviance (457,45)

AIC (465,45)

Kesimpulan

Semakin tinggi kependulian akan lingkungan alam, semakin tinggi kemungkinan memilih kendaraan listrik
Income yang lebih tinggi akan meningkatkan kemungkin memilih kendaraan listrik tetapi tidak begitu signifikan
Umur tidak memiliki pengaruh yang singnifikan

nullmod <- glm(choose_ev ~ 1, data = dataev, family = binomial)
r2_nominal <- 1- (logLik(model_nominal)/logLik(nullmod))
r2_numeric <- 1- (logLik(model_numeric_ev)/logLik(nullmod))

list(McFadden_R2_Nominal = r2_nominal,
     McFadden_R2_Numeric = r2_numeric)

## $McFadden_R2_Nominal
## 'log Lik.' 0.270279 (df=5)
## 
## $McFadden_R2_Numeric
## 'log Lik.' 0.2702601 (df=4)

Dari uji McFadden’s R Squared diatas menguji apakah model mana yang lebih baik dalam memprediksi, model nominal dan numeric tidak memiliki perbedaan nilai yang besar. Akan tetapi, model nominal memiliki angka yang sedikit lebih besar daripada model numerik.

Visualisasi Prediksi

Plot untuk Model Nominal

# Ubah environmental_awareness ke faktor (nominal)
dataev_nominal <- dataev %>%
  mutate(environmental_awareness = factor(environmental_awareness, levels = c("Low", "Medium", "High")))

# Fit model nominal
model_nominal_ev <- glm(choose_ev ~ income + age + environmental_awareness, 
                        data = dataev_nominal, 
                        family = binomial)

# Tambahkan kolom predicted
dataev_nominal <- dataev_nominal %>%
  mutate(predicted = predict(model_nominal_ev, type = "response"))

library(ggplot2)

# Plot untuk model nominal
ggplot(dataev_nominal, aes(x = income, y = predicted, color = environmental_awareness)) +
  geom_point(alpha = 0.6) +
  labs(title = "Prediksi Probabilitas Memilih Mobil Listrik\n(Environmental Awareness sebagai Nominal)",
       x = "Pendapatan (Income)",
       y = "Probabilitas Prediksi") +
  theme_minimal()

Plot untuk Model Numerik

# Tambahkan kolom predicted ke dataev_numeric
dataev_numeric <- dataev_numeric %>% mutate(predicted = predict(model_numeric_ev, type = "response"))

# Plot untuk model ordinal
ggplot(dataev_numeric, aes(x = income, y = predicted, color = as.factor(environmental_awareness_numeric))) +
  geom_point(alpha = 0.6) +
  labs(title = "Prediksi Probabilitas Memilih Mobil Listrik\n(Environmental Awareness sebagai Numeric)",
       x = "Pendapatan (Income)",
       y = "Probabilitas Prediksi",
       color = "Env Awareness\n(Numeric)") +
  theme_minimal()

Kesimpulan Akhir

Model dengan treat sebagai rasio sedikit lebih efisien untuk digunakan, karena dengan melihat AIC nya yang lebih rendah daripada AIC model dengan treat sebagai nominal

Pemilihan Model Regresi Logistik dan Evaluasi

Membangun Model Regresi Logistik : Pendekatan Confirmatori dan Exploratory

Dalam analisis regresi logistik, pemilihan model sangat krusial untuk mendapatkan model yang baik dalam memprediksi probabilitas kejadian suatu peristiwa (respon biner). Dua pendekatan utama dalam membangun model adalah pendekatan Confirmatory dan Exploratory.

Confirmatory (Pendekatan Konfirmatori)

Pendekatan ini digunakan ketika peneliti telah memiliki teori atau hipotesis yang jelas mengenai efek atau hubungan antara variabel prediktor dan respon.
- Ciri-ciri:
  - Model dibangun berdasarkan teori atau hasil penelitian sebelumnya.
  - Tujuan utamanya adalah menguji apakah efek tersebut benar-benar signifikan, bukan sekadar mencari model terbaik.
  - Peneliti biasanya menyusun model penuh terlebih dahulu, lalu menguji apakah penambahan atau pengurangan suatu efek (misalnya, interaksi) memberikan peningkatan model secara signifikan.
  - Uji signifikansi dilakukan dengan membandingkan model dengan efek tertentu dan model tanpa efek tersebut, misalnya dengan Likelihood Ratio Test.
- Contoh penggunaan: Misalnya, teori menyatakan bahwa faktor x1 dan x2 mempengaruhi probabilitas seseorang membeli produk. Maka model logistik dibangun langsung dengan x1 dan x2, lalu diuji apakah kontribusi x2 benar-benar signifikan.
Exploratory (Pendekatan Ekspolaratori)

Pendekatan ini digunakan ketika peneliti belum memiliki teori yang pasti atau ingin mengeksplorasi hubungan potensial antar variabel.
- Ciri-ciri:
  - Model dibangun secara bertahap dengan tujuan menemukan kombinasi prediktor terbaik.
  - Pemilihan variabel dilakukan berdasarkan kriteria statistik, seperti AIC, deviance, atau log-likelihood. Proses seleksi dilakukan melalui:
  - Forward Selection: Mulai dari model kosong, satu per satu variabel dimasukkan jika signifikan.
  - Backward Elimination: Mulai dari model penuh, variabel yang tidak signifikan dikeluarkan.
  - Stepwise Selection: Gabungan dari keduanya, variabel dapat masuk dan keluar secara dinamis.

Tujuan :

Menemukan model yang parsimonious, yaitu cukup sederhana namun memiliki performa prediksi yang baik.\

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

Simulasi Data

Menggunakan studi kasus “Menentukan Kemungkinan Membeli Produk Baru” menggunakan data simulasi dengan :

X1 : Sikap terhadap teknologi (Skala Kontinu)

X2 : Jenis Kelamin (0 = Perempuan, 1 = Laki-laki)

X3 : Kecenderungan mengikuti tren baru (Kontinu)

y : Apakah konsumen membeli produk baru

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

## Warning: package 'MASS' was built under R version 4.4.3

## 
## Attaching package: 'MASS'

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

 library(caret)
 library(pROC)

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

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

## 
## Attaching package: 'pROC'

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

 library(DescTools)

set.seed(64)
n <- 160
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
lin_pred <--0.5 + 1 * x1- 0.8 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
df <- data.frame(y =as.factor(y), x1, x2, x3)
head(df)

##   y          x1 x2          x3
## 1 0 -1.70796214  0 -0.94197717
## 2 0 -1.71288978  0 -0.04546181
## 3 0 -1.34170237  1  1.19686425
## 4 1  1.54077283  0  0.95864712
## 5 1  0.01518175  1  1.42538355
## 6 0 -0.32569460  0 -0.83134079

Pemilihan Model

Model Full

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

## 
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3812     0.2486  -1.533  0.12516    
## x1            1.3631     0.2663   5.119 3.07e-07 ***
## x2           -1.2831     0.4197  -3.057  0.00223 ** 
## x3            0.4769     0.2038   2.340  0.01929 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 208.39  on 159  degrees of freedom
## Residual deviance: 160.48  on 156  degrees of freedom
## AIC: 168.48
## 
## Number of Fisher Scoring iterations: 5

Metode Stepwise : Forward, Backward, dan Keduanya

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)

##               df      AIC
## model_full     4 168.4759
## step_forward   4 168.4759
## step_backward  4 168.4759
## step_both      4 168.4759

Evaluasi Model ROC dan UAC

Kurva ROC (Receiver Operating Characteristic) merupakan grafik yang menggambarkan trade off antara sensitivity (recall) dan 1 - specificity (false positive rate)

Sumbu Y : Sensitivity = True Positive Rate = TP / (TP + FN) Sumbu X : 1- Specificity = False Positive Rate = FP / (FP + TN)

Pergerakan Kurva :

Saat cut-off menurun, model mengklasifikasikan lebih banyak pengamatan sebagai positif, dalam hal ini :
- Sensitivitas naik
- Spesifisitas turun
Saat cut-off naik, model menjadi lebih konservatif :

-   Sensitivitas turun

-   Spesifisitas naik

Kurva ROC ideal memiliki bentuk seperti :

Naik tajam secara vertikal hingga mencapai sensitivitas = 1
Lalu bergerak horizontal menuju 1 - specificity = 1
Area AUC mendekati 1

Interpretasi AUC :

AUC = 0,5 : model tidak lebih baik dari tebak acak
AUC > 0,7 : model cukup baik
AUC > 0,9 : model sangat baik

Fungsi Kurva ROC :

Untuk membandingkan performa beberapa model klasifikasi
Untuk memilih threshold (cut-off) optimal berdasarkan kebutuhan aplikasi

library(pROC)
pred_prob <- predict(step_both, type = "response")
roc_obj <- roc(as.numeric(as.character(df$y)), pred_prob)

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

## Setting direction: controls < cases

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

auc(roc_obj)

## Area under the curve: 0.8126

Dari nilai AUC yang didapat yaitu 0,8126 yang dimana model ini sudah terbilang baik . Ini berarti bahwa jika kita memilih secara acak satu orang konsumen yang membeli produk baru dan satu orang yang tidak, maka model ini hanya memiliki 81,26% untuk memprediksi mana yang membeli dan mana yang tidak.

Pseudo R Squared

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

## Nagelkerke   McFadden 
##  0.3554282  0.2299385

Interpretasi Pseudo R Squared:

Nagelkerke : model menjelaskan 35,5% variabilitas atau keberagaman dalam data
McFadden : Nilai yang dihasilkan yaitu 0,23. Hasil ini menunjukkan bahwa model yang digunakan termasuk kategori yang sangat baik

Tabel Klasifikasi dan Evaluasi

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

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 88 24
##          1 15 33
##                                           
##                Accuracy : 0.7562          
##                  95% CI : (0.6822, 0.8206)
##     No Information Rate : 0.6438          
##     P-Value [Acc > NIR] : 0.001517        
##                                           
##                   Kappa : 0.4492          
##                                           
##  Mcnemar's Test P-Value : 0.200185        
##                                           
##             Sensitivity : 0.5789          
##             Specificity : 0.8544          
##          Pos Pred Value : 0.6875          
##          Neg Pred Value : 0.7857          
##              Prevalence : 0.3563          
##          Detection Rate : 0.2062          
##    Detection Prevalence : 0.3000          
##       Balanced Accuracy : 0.7167          
##                                           
##        'Positive' Class : 1               
##

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

## Sensitivity Specificity 
##   0.5789474   0.8543689

Interpretasi :

Kemampuan model mendeteksi bahwa konsumen membeli produk nya sebesar 57,9%
Kemampuan model mendeteksi bahwa konsumen tidak membeli produk nya sebesar 85,44%

Metode Perbandingan Model dalam Regresi Logistik

Pembuatan Model Bertingkat

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

Perbandingan AIC dan Deviance

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

##     Model      AIC Deviance
## 1 Model 1 180.3754 176.3754
## 2 Model 2 172.3532 166.3532
## 3 Model 3 168.4759 160.4759

Likelihood Ratio Test

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

## Analysis of Deviance Table
## 
## Model 1: y ~ x1
## Model 2: y ~ x1 + x2
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       158     176.38                        
## 2       157     166.35  1   10.022 0.001547 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

## Analysis of Deviance Table
## 
## Model 1: y ~ x1 + x2
## Model 2: y ~ x1 + x2 + x3
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       157     166.35                       
## 2       156     160.48  1   5.8773  0.01534 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prinsip Parsimony

Memang, model yang kompleks sering memiliki AIC dan deviance yang lebih kecil. Namun dalam prinsip parsimony memerhatikan hal hal berikut :

Model sederhana agar lebih mudah diinterpretasikan.
Jika penurunan AIC tidak signifikan, pilih model lebih sederhana.
Parsimony mencegah overfitting.

Precision-Recall Curve (PR Curve)

Kurva Precision-Recall (PR) adalah alat evaluasi performa model klasifikasi, khususnya sangat berguna saat bekerja dengan data yang tidak seimbang (class imbalance).

Regresi Logistik Multinomial

Distribusi Multinomial

Distribusi multinomial adalah perluasan dari distribusi binomial untuk lebih dari dua kategori. Jika X1,X2,…,X𝑘 menyatakan banyaknya kejadian dalam masing-masing dari 𝑘 kategori, maka:

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

dengan \[ \sum_{i=1}^{k} x_i = n \quad \text{dan} \quad \sum_{i=1}^{k} p_i = 1. \]

Regresi Logistik Multinomial

Regresi Logistik Multinomial adalah perluasan dari regresi logistik. Model Regresi Logistik Multinomial digunaka untuk memodelkan hubungan antara satu variabel respon kategorik (>2 kategori) dan satu atau leih variabel prediktor.

\[ \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,…, K-1

Estimasi Parameter

Estimasi dilakukan menggunakan metode maximum likelihood dengan algoritma iteratif seperti Newton Raphson.

Log-Likelihood :

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

\[ \text{dengan } \pi_{ij} = P(Y_i = j \mid \mathbf{x}_i) \text{ dan } y_{ij} = 1 \text{ jika } Y_i = j. \]

Contoh Kasus : Simulasi Data

Judul : Preferensi Transportasi Karyawan ke Kantor

Latar Belakang:

Sebuah perusahaan ingin memahami faktor-faktor yang memengaruhi pilihan moda transportasi yang digunakan karyawan untuk pergi ke kantor. Perusahaan melakukan survei terhadap 200 karyawan dan mengumpulkan data berikut:

Moda Transportasi : Mobil Pribadi, Sepeda Motor, Transportasi Umum
Usia
Income Level : Rendah, Menengah, Tinggi
Jarak rumah ke Kantor (dalam km)

Tujuan Analisis :

Mengetahui bagaimana usia, tingkat pendapatan, dan jarak rumah ke kantor memengaruhi pilihan moda transportasi

Simulasi Data

library(nnet)
set.seed(164)
n <- 200

# Variabel independen
IncomeLevel <- sample(c("Rendah", "Menengah", "Tinggi"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 35, sd = 7))
Distance <- round(pmax(rnorm(n, mean = 10, sd = 5), 0.5), 1)  # Jarak dalam km

# Simulasikan TransportMode berdasarkan probabilitas berbeda per IncomeLevel
TransportMode <- sapply(IncomeLevel, function(inc) {
  if (inc == "Rendah") {
    sample(c("Transportasi Umum", "Sepeda Motor", "Mobil Pribadi"), size = 1,
           prob = c(0.6, 0.3, 0.1))
  } else if (inc == "Menengah") {
    sample(c("Transportasi Umum", "Sepeda Motor", "Mobil Pribadi"), size = 1,
           prob = c(0.3, 0.4, 0.3))
  } else {
    sample(c("Transportasi Umum", "Sepeda Motor", "Mobil Pribadi"), size = 1,
           prob = c(0.1, 0.3, 0.6))
  }
})

# Bentuk data frame
data_reglogmul <- data.frame(
  TransportMode = factor(TransportMode),
  Age,
  IncomeLevel = factor(IncomeLevel),
  Distance
)

# Set "Transportasi Umum" sebagai baseline
data_reglogmul$TransportMode <- relevel(data_reglogmul$TransportMode, ref = "Transportasi Umum")

# Tampilkan sebagian data
head(data_reglogmul)

##       TransportMode Age IncomeLevel Distance
## 1     Mobil Pribadi  31    Menengah      7.5
## 2 Transportasi Umum  32    Menengah     14.9
## 3     Mobil Pribadi  32      Tinggi     11.8
## 4 Transportasi Umum  45    Menengah      9.1
## 5     Mobil Pribadi  22      Rendah     14.2
## 6     Mobil Pribadi  33      Tinggi     13.9

Estimasi Model

model_mnlogit <- multinom(TransportMode ~ IncomeLevel + Age + Distance, data = data_reglogmul)

## # weights:  18 (10 variable)
## initial  value 219.722458 
## iter  10 value 196.944744
## final  value 196.649721 
## converged

summary(model_mnlogit)

## Call:
## multinom(formula = TransportMode ~ IncomeLevel + Age + Distance, 
##     data = data_reglogmul)
## 
## Coefficients:
##               (Intercept) IncomeLevelRendah IncomeLevelTinggi         Age
## Mobil Pribadi    1.608988         -2.417500         0.6421718 -0.04385405
## Sepeda Motor     1.670733         -1.223187         0.2537100 -0.04149782
##                 Distance
## Mobil Pribadi 0.03817385
## Sepeda Motor  0.02639670
## 
## Std. Errors:
##               (Intercept) IncomeLevelRendah IncomeLevelTinggi        Age
## Mobil Pribadi    1.133280         0.5711007         0.4784898 0.02870238
## Sepeda Motor     1.046397         0.4494412         0.4882161 0.02652377
##                 Distance
## Mobil Pribadi 0.04400823
## Sepeda Motor  0.03972819
## 
## Residual Deviance: 393.2994 
## AIC: 413.2994

Nilai P-Value dan Interpretasi

z <- summary(model_mnlogit)$coefficients / summary(model_mnlogit)$standard.errors
pval <- 2 * (1- pnorm(abs(z)))
round(pval, 4)

##               (Intercept) IncomeLevelRendah IncomeLevelTinggi    Age Distance
## Mobil Pribadi      0.1557            0.0000            0.1796 0.1265   0.3857
## Sepeda Motor       0.1103            0.0065            0.6033 0.1177   0.5064

Interpretasi :

Model ini membandingkan dua kateogri terhadap baseline yaitu Transportasi Umum
Untuk mobil pribadi dibandingkan dengan transportasi umum, variabel-variabel yang ada menunjukkan p-value > 0,05 kecuali variabel IncomeLevelRendah yang memiliki p-value = 0,0000 yang dimana <0,05
Untuk sepeda motor dibandingkan dengan transportasi umum, variabel-variabel yang ada menunjukkan p-value > 0,05 kecuali variabel IncomeLevelRendah yang memiliki p-value = 0,0065 yang dimana <0,05
IncomeLevelRendah merupakan satu-satunya variabel yang konsisten signifikan sehingga dapat disimpulkan bahwa karyawan yang memiliki penghasilan rendah cenderung lebih mungkin memilih sepeda motor dan lebih kecil kemungkinan memilih mobil pribadi dibanding karyawan yang berpenghasilan menengah

Prediksi dan Validasi

data_reglogmul$Predicted <- predict(model_mnlogit)
table(Predicted = data_reglogmul$Predicted, Actual = data_reglogmul$TransportMode)

##                    Actual
## Predicted           Transportasi Umum Mobil Pribadi Sepeda Motor
##   Transportasi Umum                40             8           23
##   Mobil Pribadi                    17            39           29
##   Sepeda Motor                     11            17           16

# 1. Buat kombinasi semua level prediktor (misalnya kamu punya 3 variabel)
new_data <- expand.grid(
  Age = mean(data_reglogmul$Age),  # bisa juga gunakan range jika ingin
  IncomeLevel = levels(data_reglogmul$IncomeLevel),
  Distance = mean(data_reglogmul$Distance)
)

# 2. Prediksi probabilitas dari model multinomial
predicted_probs <- predict(model_mnlogit, newdata = new_data, type = "probs")

# 3. Ubah nama kolom probabilitas agar jelas
colnames(predicted_probs) <- paste0("prob.", colnames(predicted_probs))

# 4. Gabungkan dengan data prediktor
result_df <- cbind(new_data, predicted_probs)

# 5. (Opsional) Tampilkan
print(result_df)

##      Age IncomeLevel Distance prob.Transportasi Umum prob.Mobil Pribadi
## 1 34.955    Menengah   9.7785              0.2391989         0.37489568
## 2 34.955      Rendah   9.7785              0.6193861         0.08653782
## 3 34.955      Tinggi   9.7785              0.1650694         0.49171037
##   prob.Sepeda Motor
## 1         0.3859054
## 2         0.2940761
## 3         0.3432203

Regresi Logistik Ordinal

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

Model ini berbeda dengan:

Regresi logistik biner : hanya terdapat 2 kategori
Regresi logistik multinomial : memiliki kategori >2 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 :

𝛼𝑗: intercept khusus untuk kategori ke-j
𝛽: koefisien regresi (sama untuk semua kategori kumulatif)

Untuk 𝑐 kategori, terdapat (𝑐 − 1) model logit kumulatif.

Interpretasi Koefisien

Koefisien 𝛽 menjelaskan efek 𝑥 terhadap kemungkinan berada pada kategori yang lebih rendah atau sama. Jika 𝛽>0: semakin besar 𝑥, semakin tinggi peluang berada di kategori rendah. Jika 𝛽<0: semakin besar 𝑥, semakin besar peluang berada di kategori tinggi.

Odds ratio:

\[ \text{OR} = e^{\beta} \]

Contoh Kasus : Simulasi Data

Contoh Kasus : Kepuasan terhadap Layanan Internet

library(tibble)
library(dplyr)
library(tidyr)

# Data ordinal buatan
internet_satisfaction <- tribble(
  ~Gender, ~AgeGroup, ~Provider, ~Satisfaction, ~Count,
  "Male",   "18-25", "ISP_A", "Very Dissatisfied", 10,
  "Male",   "18-25", "ISP_A", "Dissatisfied",      20,
  "Male",   "18-25", "ISP_A", "Neutral",           40,
  "Male",   "18-25", "ISP_A", "Satisfied",         60,
  "Male",   "18-25", "ISP_A", "Very Satisfied",    30,
  
  "Female", "26-35", "ISP_B", "Very Dissatisfied", 5,
  "Female", "26-35", "ISP_B", "Dissatisfied",      15,
  "Female", "26-35", "ISP_B", "Neutral",           30,
  "Female", "26-35", "ISP_B", "Satisfied",         50,
  "Female", "26-35", "ISP_B", "Very Satisfied",    45,
  
  "Male",   "36-50", "ISP_B", "Very Dissatisfied", 2,
  "Male",   "36-50", "ISP_B", "Dissatisfied",      5,
  "Male",   "36-50", "ISP_B", "Neutral",           20,
  "Male",   "36-50", "ISP_B", "Satisfied",         40,
  "Male",   "36-50", "ISP_B", "Very Satisfied",    60,
  
  "Female", "18-25", "ISP_A", "Very Dissatisfied", 7,
  "Female", "18-25", "ISP_A", "Dissatisfied",      10,
  "Female", "18-25", "ISP_A", "Neutral",           35,
  "Female", "18-25", "ISP_A", "Satisfied",         55,
  "Female", "18-25", "ISP_A", "Very Satisfied",    40
)

# Ubah ke tipe faktor ordinal & kategorik
internet_satisfaction <- internet_satisfaction %>%
  mutate(
    Satisfaction = factor(Satisfaction, 
                          levels = c("Very Dissatisfied", "Dissatisfied", 
                                     "Neutral", "Satisfied", "Very Satisfied"),
                          ordered = TRUE),
    Gender = factor(Gender),
    AgeGroup = factor(AgeGroup),
    Provider = factor(Provider)
  )

Ekspansi Data

expanded_data <- uncount(internet_satisfaction, weights = Count)
head(expanded_data)

## # A tibble: 6 × 4
##   Gender AgeGroup Provider Satisfaction     
##   <fct>  <fct>    <fct>    <ord>            
## 1 Male   18-25    ISP_A    Very Dissatisfied
## 2 Male   18-25    ISP_A    Very Dissatisfied
## 3 Male   18-25    ISP_A    Very Dissatisfied
## 4 Male   18-25    ISP_A    Very Dissatisfied
## 5 Male   18-25    ISP_A    Very Dissatisfied
## 6 Male   18-25    ISP_A    Very Dissatisfied

Estimasi Model Ordinal

model_ordinal <- polr(Satisfaction ~ Gender + AgeGroup + Provider, data = expanded_data, Hess = TRUE)

## Warning in polr(Satisfaction ~ Gender + AgeGroup + Provider, data =
## expanded_data, : design appears to be rank-deficient, so dropping some coefs

summary(model_ordinal)

## Call:
## polr(formula = Satisfaction ~ Gender + AgeGroup + Provider, data = expanded_data, 
##     Hess = TRUE)
## 
## Coefficients:
##                  Value Std. Error t value
## GenderMale    -0.41360     0.2053 -2.0146
## AgeGroup26-35  0.08265     0.2125  0.3889
## AgeGroup36-50  1.22257     0.2224  5.4982
## 
## Intercepts:
##                                Value    Std. Error t value 
## Very Dissatisfied|Dissatisfied  -3.1344   0.2458   -12.7532
## Dissatisfied|Neutral            -1.9004   0.1792   -10.6064
## Neutral|Satisfied               -0.5970   0.1546    -3.8617
## Satisfied|Very Satisfied         0.9490   0.1584     5.9895
## 
## Residual Deviance: 1594.057 
## AIC: 1608.057

Uji Signifikansi dan P-Value

coefs <- coef(summary(model_ordinal))
p_values <- pnorm(abs(coefs[, "t value"]), lower.tail = FALSE) * 2
cbind(coefs, "p value" = round(p_values, 4)) %>%
  knitr::kable(caption = "Koefisien dan p-value dari Model Regresi Logistik Ordinal")

Koefisien dan p-value dari Model Regresi Logistik Ordinal
	Value	Std. Error	t value	p value
GenderMale	-0.4136006	0.2053000	-2.0146162	0.0439
AgeGroup26-35	0.0826455	0.2124930	0.3889329	0.6973
AgeGroup36-50	1.2225706	0.2223569	5.4982347	0.0000
Very Dissatisfied\|Dissatisfied	-3.1344008	0.2457735	-12.7532073	0.0000
Dissatisfied\|Neutral	-1.9004174	0.1791766	-10.6063951	0.0000
Neutral\|Satisfied	-0.5970351	0.1546045	-3.8616928	0.0001
Satisfied\|Very Satisfied	0.9489725	0.1584393	5.9895024	0.0000

Kesimpulan

Koefisien :
- Semakin besar koefisien nya, maka akan cenderung puas akan layanan internet. Begitupun sebaliknya.
- GenderMale : DIkarenakan p-value = 0,0436 yang mana <0,5 sehingga dapatdisimpulkan bahwa laki-laki cenderung memiliki kepuasan lebih rendah dibanding perempuan
- AgeGroup26-35 : Dikarenakan p-value = 0,6973 yang mana >0,5, kepuasan kelompok usia ini tidak jauh berbeda dengan kelompok usia 18-25 tahun
- AgeGroup36-50 : Dikarenakan p-value = 0,0001 yang mana <0,5 dan juga melihat koefisien nya yaitu 1,2226. Dapat disimpulkan bahwa kelompok usia ini cenderung sangat puas akan layan internet yang didapatkan.

Log Linear Model

Analisis data kategorik merupakan bagian penting dalam statistika terapan karena banyak fenomena di dunia nyata yang menghasilkan data dalam bentuk kategori, seperti jenis kelamin, status pekerjaan, tingkat pendidikan, preferensi konsumen, atau hasil diagnosis medis. Data kategori ini umumnya dianalisis menggunakan tabel kontingensi, model log-linier, dan model regresi logistik. Masing-masing pendekatan memiliki kekuatan dan kelemahan tergantung pada tujuan analisis dan struktur data.

Tabel kontingensi digunakan sebagai langkah awal eksplorasi untuk melihat hubungan antara dua atau lebih variabel kategorik. Misalnya, dalam studi tentang efek obat terhadap serangan jantung, tabel kontingensi dapat menyajikan jumlah pasien yang mengalami atau tidak mengalami serangan jantung, berdasarkan jenis obat yang dikonsumsi. Tabel ini membantu mengidentifikasi pola awal dan menghitung ukuran asosiasi seperti odds ratio, risk ratio, dan chi-square statistic untuk menguji independensi antar variabel.

Namun, ketika ingin membangun model statistik yang dapat mengendalikan efek dari banyak variabel dan interaksinya secara simultan, maka model log-linier menjadi sangat berguna. Model log-linier adalah bentuk khusus dari Generalized Linear Model (GLM) yang digunakan pada frekuensi sel dalam tabel kontingensi dan mengasumsikan distribusi Poisson. Tidak seperti regresi logistik, model log-linier tidak menetapkan variabel mana yang dependen dan mana yang independen, karena seluruh variabel diperlakukan secara simetris. Model ini lebih cocok ketika tujuan analisis adalah untuk memahami struktur asosiasi atau independensi antar variabel, bukan untuk prediksi.

Struktur model log-linier ditentukan berdasarkan efek utama dari masing-masing variabel serta interaksi di antara variabel-variabel tersebut. Misalnya, dalam tabel kontingensi tiga arah (misalnya: jenis kelamin, status merokok, dan penyakit paru), model log-linier dapat membedakan apakah interaksi dua variabel cukup menjelaskan data, ataukah diperlukan interaksi tiga arah untuk menangkap struktur asosiasinya. Penyesuaian model dapat dilakukan menggunakan metode likelihood ratio test untuk membandingkan model sederhana dengan model lebih kompleks.

Di sisi lain, regresi logistik adalah pendekatan paling umum ketika terdapat satu variabel kategorik yang secara eksplisit dianggap sebagai variabel dependen (misalnya, kejadian penyakit: ya/tidak), dan satu atau lebih variabel kategorik atau numerik sebagai prediktor. Model ini memodelkan logit dari probabilitas kejadian (yaitu log odds), dan sangat berguna dalam studi observasional dan eksperimental untuk menjelaskan atau memprediksi peluang suatu outcome. Regresi logistik juga memiliki ekstensi untuk outcome kategorik lebih dari dua kelas, seperti regresi logistik multinomial dan regresi logistik ordinal.

Dengan demikian, meskipun ketiganya beroperasi pada data kategorik, tabel kontingensi bersifat deskriptif, model log-linier bersifat eksploratif terhadap hubungan simetris, sedangkan regresi logistik bersifat prediktif terhadap outcome kategorik. Pemilihan metode tergantung pada apakah fokus utama analisis adalah deskripsi, eksplorasi struktur, atau prediksi hasil berdasarkan variabel penjelas. Kombinasi dari ketiganya sering digunakan dalam praktik untuk memperoleh pemahaman komprehensif dari data kategorik yang dianalisis

Ringkasan dalam analisis data kategori tedapat beberapa pendekatan statistik yang umum digunakan antara lain :

Tabel Kontingensi: penyajian frekuensi gabungan dari dua atau lebih variabel kategorik.
Model Loglinear: digunakan untuk memodelkan struktur asosiasi di dalam tabel kontingensi tanpa menganggap ada variabel dependen.
Model Regresi Logistik: digunakan untuk memodelkan probabilitas dari kategori variabel dependen berdasarkan variabel independen.

Tabel Kontingensi

Tabel kontingensi menyajikan frekuensi dari kombinasi antar variabel. Misal :

data2x2<- matrix(c(189, 10845, 104, 10933), nrow = 2, byrow = TRUE)
dimnames(data2x2) <- list("Kelompok" = c("Placebo", "Aspirin"), "Kejadian" = c("Serangan Jantung (Ya)", "Serangan Jantung (Tidak)"))
print(data2x2)

##          Kejadian
## Kelompok  Serangan Jantung (Ya) Serangan Jantung (Tidak)
##   Placebo                   189                    10845
##   Aspirin                   104                    10933

Model Log Linear

Model loglinear memodelkan logaritma dari ekspektasi frekuensi sel dalam tabel kontingensi.

\[ \log(\mu_{ij}) = \mu + \lambda_i^A + \lambda_j^B + \lambda_{ij}^{AB} \]

Cocok untuk menyelidiki asosiasi dan independensi antar variabel.
Tidak membedakan antara variabel respon dan penjelas.
Umumnya digunakan pada tabel multi-dimensi (2 arah, 3 arah, dst).

library(MASS)
loglm(~ Kelompok * Kejadian, data = data2x2)

## Call:
## loglm(formula = ~Kelompok * Kejadian, data = data2x2)
## 
## Statistics:
##                  X^2 df P(> X^2)
## Likelihood Ratio   0  0        1
## Pearson            0  0        1

Model Regresi Logistik

Model Regresi Logistik Biner :

\[ \log\left(\frac{p}{1 - p}\right) = \beta_0 + \beta_1 x \]

Digunakan jika ada variabel dependen kategorik (biasanya biner).
Bertujuan untuk memprediksi probabilitas suatu outcome.
Umumnya digunakan dalam studi observasional atau eksperimental.

Contoh :

data_glm <- data.frame(
Kejadian = c(1, 0, 1, 0),
Kelopok = factor(c("Placebo", "Placebo", "Aspirin", "Aspirin")),
Frek = c(189, 10845, 104, 10933))

model_logit <- glm(Kejadian ~ Kelopok, weights = Frek, family = binomial, data = data_glm)
summary(model_logit)

## 
## Call:
## glm(formula = Kejadian ~ Kelopok, family = binomial, data = data_glm, 
##     weights = Frek)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -4.65515    0.09852 -47.250  < 2e-16 ***
## KelopokPlacebo  0.60544    0.12284   4.929 8.28e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3114.7  on 3  degrees of freedom
## Residual deviance: 3089.3  on 2  degrees of freedom
## AIC: 3093.3
## 
## Number of Fisher Scoring iterations: 7

Tabel Kontingensi dan Model Loglinier

Tabel kontingensi menyajikan frekuensi dari kombinasi kategori antar dua atau lebih variabel. Misal :

matrix(c(189, 104,10845, 10933), nrow=2,
 dimnames = list(Kelompok = c("Placebo", "Aspirin"),
 Kejadian = c("Ya", "Tidak")))

##          Kejadian
## Kelompok   Ya Tidak
##   Placebo 189 10845
##   Aspirin 104 10933

Model log-linier untuk tabel I x J dapat dituliskan sebagai :

\[ \log(\mu_{ij}) = \mu + \lambda_i^T + \lambda_j^R + \lambda_{ij}^{TR} \]

Model Saturated

Model saturated atau model penuh menyertakan seluruh efek utama dan interaksi. Model ini biasanya digunakan sebagai pembanding saat menguji model yang lebih sederhana misalnya menggunakan Likelihood Ratio Test.

library(MASS)
data2x <- matrix(c(189, 10845, 104, 10933), nrow=2, byrow=TRUE)
dimnames(data2x) <- list(Kelompok = c("Placebo", "Aspirin"), Kejadian = c("Ya", "Tidak"))
ftable(data2x)

##          Kejadian    Ya Tidak
## Kelompok                     
## Placebo             189 10845
## Aspirin             104 10933

Model saturated dapat dipasang dengan loglm dari package {MASS}:

model_saturated <- loglm(~ Kelompok * Kejadian, data = data2x)
summary(model_saturated)

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

Model Independent

Model independen mengasumsikan bahwa tidak ada interaksi antara variabel :

\[ \log(\mu_{ij}) = \mu + \lambda_i^T + \lambda_j^R \]

Model ini menguji hipotesis bahwa variabel X dan Y saling independen.

model_indep <- loglm(~ Kelompok + Kejadian, data = data2x)
summary(model_indep)

## Formula:
## ~Kelompok + Kejadian
## attr(,"variables")
## list(Kelompok, Kejadian)
## attr(,"factors")
##          Kelompok Kejadian
## Kelompok        1        0
## Kejadian        0        1
## attr(,"term.labels")
## [1] "Kelompok" "Kejadian"
## 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.37196  1 4.727402e-07
## Pearson          25.01388  1 5.691897e-07

Estimasi Parameter

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

## [1] 0.6054377

Model Lebih Sederhana dan Perbandingan Model

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

anova(model_indep, model_saturated)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Kelompok + Kejadian 
## Model 2:
##  ~Kelompok * Kejadian 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   25.37196  1                                    
## Model 2    0.00000  0   25.37196         1              0
## Saturated  0.00000  0    0.00000         0              1

Studi Kasus : Simulasi Data

Judul : Hubungan antara Status Merokok dan Status Kesehatan

data_merokok <- matrix(c(60,180,110,
                         30,240,220),
                       nrow = 2, byrow = TRUE,
                       dimnames = list(
                         Merokok = c("Perokok", "Bukan Perokok"),
                         Kesehatan = c("Buruk", "Cukup", "Baik")
                       ))

ftable(data_merokok)

##               Kesehatan Buruk Cukup Baik
## Merokok                                 
## Perokok                    60   180  110
## Bukan Perokok              30   240  220

loglm(~ Merokok + Kesehatan, data = data_merokok)

## Call:
## loglm(formula = ~Merokok + Kesehatan, data = data_merokok)
## 
## Statistics:
##                       X^2 df     P(> X^2)
## Likelihood Ratio 32.72998  2 7.812242e-08
## Pearson          32.81633  2 7.482130e-08

Interpretasi :

Dikarenakan p-value < 0,05 maka terdapat hubungan antar variabel yang menyatakan bahwa kelompok merokok dan kesehatan itu tidak independen

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 variabel kategorik yang disajikan dalam tabel kontingensi. Model ini mengasumsikan bahwa logaritma dari nilai ekspektasi frekuensi sel (𝜇𝑖𝑗) 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)

Bentuk Log Linear

Model log-linear pada tabel 2x2 :

\[ \log(\mu_{ij}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_{ij}^{AB} \]

dengan constraint sum-to-zero :

\[ \sum_i \lambda_i^A = 0, \quad \sum_j \lambda_j^B = 0, \quad \sum_{i,j} \lambda_{ij}^{AB} = 0 \]

Analisis Data Tabel Kontingensi 2x2 : Contoh Kasus

matrix(c(189, 104,10845, 10933), nrow=2,
dimnames = list(Kelompok = c("Placebo", "Aspirin"),
Kejadian = c("Ya", "Tidak")))

##          Kejadian
## Kelompok   Ya Tidak
##   Placebo 189 10845
##   Aspirin 104 10933

Hitung Odds Ratio dan Interval Kepercayaan

# Angka-angka dalam tabel
n11 <- 189
n12 <- 10845
n21 <- 104
n22 <- 10933

# Hitung Odds Ratio
OR <- (n11 * n22) / (n12 * n21)
OR

## [1] 1.832054

# Log Odds Ratio
log_OR <- log(OR)
log_OR

## [1] 0.6054377

# Hitung Standard Error (SE)
SE <- sqrt(1/n11 + 1/n12 + 1/n21 + 1/n22)
SE

## [1] 0.1228416

# Hitung Confidence Interval untuk log(OR)
z <- 1.96  # Z-score untuk 95% CI
lower_log <- log_OR - z * SE
upper_log <- log_OR + z * SE

# Transformasi balik untuk mendapatkan CI dari OR
lower_OR <- exp(lower_log)
upper_OR <- exp(upper_log)

# Tampilkan hasil akhir
cat("OR =", round(OR, 2), "\n")

## OR = 1.83

cat("95% Confidence Interval: (", round(lower_OR, 2), ",", round(upper_OR, 2), ")\n")

## 95% Confidence Interval: ( 1.44 , 2.33 )

# Hitung Odds Ratio
OR <- (n11 * n22) / (n12 * n21)
OR

## [1] 1.832054

# Log Odds Ratio
log_OR <- log(OR)
log_OR

## [1] 0.6054377

# Hitung Standard Error (SE)
SE <- sqrt(1/n11 + 1/n12 + 1/n21 + 1/n22)
SE

## [1] 0.1228416

# Hitung Confidence Interval untuk log(OR)
z <- 1.96  # Z-score untuk 95% CI
lower_log <- log_OR - z * SE
upper_log <- log_OR + z * SE

# Transformasi balik untuk mendapatkan CI dari OR
lower_OR <- exp(lower_log)
upper_OR <- exp(upper_log)

# Tampilkan hasil akhir
cat("OR =", round(OR, 2), "\n")

## OR = 1.83

cat("95% Confidence Interval: (", round(lower_OR, 2), ",", round(upper_OR, 2), ")\n")

## 95% Confidence Interval: ( 1.44 , 2.33 )

Fitting Model Log-Linear

library(MASS)
data2x <- matrix(c(189, 10845, 104, 10933), nrow=2, byrow=TRUE)
dimnames(data2x) <- list(Kelompok = c("Placebo", "Aspirin"), Kejadian = c("Ya", "Tidak"))
ftable(data2x)

##          Kejadian    Ya Tidak
## Kelompok                     
## Placebo             189 10845
## Aspirin             104 10933

data2kali2 <-as.data.frame(as.table(data2x))
colnames(data2kali2) <-c("Kelompok", "Kejadian_Serangan_Jantung", "Freq")
data2kali2

##   Kelompok Kejadian_Serangan_Jantung  Freq
## 1  Placebo                        Ya   189
## 2  Aspirin                        Ya   104
## 3  Placebo                     Tidak 10845
## 4  Aspirin                     Tidak 10933

# Model tanpa interaksi
fit_no_inter <- glm(Freq ~ Kelompok + Kejadian_Serangan_Jantung , family = poisson, data = data2kali2)
summary(fit_no_inter)

## 
## Call:
## glm(formula = Freq ~ Kelompok + Kejadian_Serangan_Jantung, family = poisson, 
##     data = data2kali2)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    4.9868895  0.0588072   84.80   <2e-16 ***
## KelompokAspirin                0.0002718  0.0134623    0.02    0.984    
## Kejadian_Serangan_JantungTidak 4.3084830  0.0588123   73.26   <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: 27507.580  on 3  degrees of freedom
## Residual deviance:    25.372  on 1  degrees of freedom
## AIC: 67.203
## 
## Number of Fisher Scoring iterations: 4

# Model dengan interaksi
fit_inter <- glm(Freq ~ Kelompok * Kejadian_Serangan_Jantung, family = poisson, data = data2kali2)
summary(fit_inter)

## 
## Call:
## glm(formula = Freq ~ Kelompok * Kejadian_Serangan_Jantung, family = poisson, 
##     data = data2kali2)
## 
## Coefficients:
##                                                Estimate Std. Error z value
## (Intercept)                                     5.24175    0.07274  72.062
## KelompokAspirin                                -0.59736    0.12209  -4.893
## Kejadian_Serangan_JantungTidak                  4.04971    0.07337  55.195
## KelompokAspirin:Kejadian_Serangan_JantungTidak  0.60544    0.12284   4.929
##                                                Pr(>|z|)    
## (Intercept)                                     < 2e-16 ***
## KelompokAspirin                                9.95e-07 ***
## Kejadian_Serangan_JantungTidak                  < 2e-16 ***
## KelompokAspirin:Kejadian_Serangan_JantungTidak 8.28e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  2.7508e+04  on 3  degrees of freedom
## Residual deviance: -3.2951e-13  on 0  degrees of freedom
## AIC: 43.831
## 
## Number of Fisher Scoring iterations: 2

Interpretasi :

Intercept menunjukkan rata-rata log frekuensi sel yaitu 5,242
Interaksi signifikan yang mana menunjukkan bahwa adanya asosiasi antara kelompok dan kejadian serangan jantung

Analisis Data tabel Kontingensi 2x3 : Contoh Kasus

data_merokok <- matrix(c(60,180,110,
                         30,240,220),
                       nrow = 2, byrow = TRUE,
                       dimnames = list(
                         Merokok = c("Perokok", "Bukan Perokok"),
                         Kesehatan = c("Buruk", "Cukup", "Baik")
                       ))

ftable(data_merokok)

##               Kesehatan Buruk Cukup Baik
## Merokok                                 
## Perokok                    60   180  110
## Bukan Perokok              30   240  220

# Ubah menjadi data.frame untuk glm
data2x3 <- as.data.frame(as.table(data_merokok))
colnames(data2x3) <- c("Meroko", "Kesehatan", "Freq")
data2x3

##          Meroko Kesehatan Freq
## 1       Perokok     Buruk   60
## 2 Bukan Perokok     Buruk   30
## 3       Perokok     Cukup  180
## 4 Bukan Perokok     Cukup  240
## 5       Perokok      Baik  110
## 6 Bukan Perokok      Baik  220

Fitting Model Log Linear

# Model log-linear tanpa interaksi (asumsi independen)
fit_no_inter <- glm(Freq ~ Meroko + Kesehatan, family = poisson, data = data2x3)
summary(fit_no_inter)

## 
## Call:
## glm(formula = Freq ~ Meroko + Kesehatan, family = poisson, data = data2x3)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          3.62434    0.11304  32.063  < 2e-16 ***
## MerokoBukan Perokok  0.33647    0.06999   4.808 1.53e-06 ***
## KesehatanCukup       1.54045    0.11615  13.262  < 2e-16 ***
## KesehatanBaik        1.29928    0.11892  10.926  < 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: 300.91  on 5  degrees of freedom
## Residual deviance:  32.73  on 2  degrees of freedom
## AIC: 80.033
## 
## Number of Fisher Scoring iterations: 4

# Model log-linear dengan interaksi (untuk cek asosiasi)
fit_inter <- glm(Freq ~ Meroko * Kesehatan, family = poisson, data = data2x3)
summary(fit_inter)

## 
## Call:
## glm(formula = Freq ~ Meroko * Kesehatan, family = poisson, data = data2x3)
## 
## Coefficients:
##                                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          4.0943     0.1291  31.715  < 2e-16 ***
## MerokoBukan Perokok                 -0.6931     0.2236  -3.100 0.001936 ** 
## KesehatanCukup                       1.0986     0.1491   7.370 1.71e-13 ***
## KesehatanBaik                        0.6061     0.1605   3.777 0.000159 ***
## MerokoBukan Perokok:KesehatanCukup   0.9808     0.2444   4.014 5.98e-05 ***
## MerokoBukan Perokok:KesehatanBaik    1.3863     0.2523   5.495 3.90e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3.0091e+02  on 5  degrees of freedom
## Residual deviance: 8.1712e-14  on 0  degrees of freedom
## AIC: 51.303
## 
## Number of Fisher Scoring iterations: 3

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

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

Model ini memuat semua kemungkinan interaksi, termasuk interaksi tiga arah (X, Y, dan Z).

Model Homogen

\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} + \lambda_{ik}^{XZ} + \lambda_{jk}^{YZ} \]
Model ini hanya mengakomodasi interaksi dua arah antar variabel tanpa memasukkan interaksi tiga arah.

Model Conditional

Conditional pada X

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

Memuat interaksi X dengan Y dan X dengan Z.
Conditional pada Y

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

Memuat interaksi Y dengan X dan Y dengan Z.
Conditional pada Z

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

Memuat interaksi Z dengan X dan Z dengan Y.

Model Joint Independence

Independensi antar X & Y

\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ij}^{XY} \]
Independensi antar X & Z

\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z + \lambda_{ik}^{XZ} \]
Independensi antar Y & Z

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

Model tanpa Interaksi

\[ \log(\mu_{ijk}) = \lambda + \lambda_i^X + \lambda_j^Y + \lambda_k^Z \]

Model ini hanya memasukkan efek utama tanpa interaksi antar variabel.

Pengujian Interaksi dalam Model Log-Linear Tiga Arah

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

Pengujian Interaksi Tiga Arah (XYZ)
- Membandingkan model saturated dengan model homogenous
Pengujian Interaksi Dua Arah (XY, XZ, YZ)
- Membandingkan model homogenous dengan model conditional
- Membandingkan model conditional dengan model joint independence
- Membandingkan model join 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.

Analisis Log-Linear untuk Tabel Tiga Arah

Judul :

Analisis Sikap Masyarakat terhadap Penggunaan Transportasi Umum Berdasarkan Jenis Kelamin dan Pendidikan

Data :

# Paket yang digunakan
library(epitools)
library(DescTools)
library(lawstat)

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

# Variabel faktor
z.edu <- factor(rep(c("SMA", "Sarjana", "Pascasarjana"), each = 4))
x.sex <- factor(rep(c("Laki-laki", "Perempuan"), each = 2, times = 3))
y.att <- factor(rep(c("Yes", "No"), times = 6))

# Data fiktif: Frekuensi
counts <- c(
  # SMA
  85, 75,     # Laki-laki
  78, 82,     # Perempuan
  # Sarjana
  140, 60,    # Laki-laki
  135, 55,    # Perempuan
  # Pascasarjana
  120, 30,    # Laki-laki
  125, 25     # Perempuan
)

# Buat data frame
data_loglinear <- data.frame(
  Pendidikan     = z.edu,
  Jenis_Kelamin  = x.sex,
  Sikap          = y.att,
  Frekuensi      = counts
)

# Tampilkan data
data_loglinear

##      Pendidikan Jenis_Kelamin Sikap Frekuensi
## 1           SMA     Laki-laki   Yes        85
## 2           SMA     Laki-laki    No        75
## 3           SMA     Perempuan   Yes        78
## 4           SMA     Perempuan    No        82
## 5       Sarjana     Laki-laki   Yes       140
## 6       Sarjana     Laki-laki    No        60
## 7       Sarjana     Perempuan   Yes       135
## 8       Sarjana     Perempuan    No        55
## 9  Pascasarjana     Laki-laki   Yes       120
## 10 Pascasarjana     Laki-laki    No        30
## 11 Pascasarjana     Perempuan   Yes       125
## 12 Pascasarjana     Perempuan    No        25

Uji Model Interaksi Tiga Arah (Saturated vs Homogenous)

Penentuan Kategori Referensi

#Penentuan Kategori Reference

x.sex <- relevel(x.sex, ref = "Perempuan")
y.att <- relevel(y.att, ref = "No")
z.edu <- relevel(z.edu, ref = "Pascasarjana")

Model Saturated Model log-linear saturated memasukkan semua interkasi hingga tiga arah :

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

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

## 
## Call:
## glm(formula = counts ~ x.sex + y.att + z.edu + x.sex * y.att + 
##     x.sex * z.edu + y.att * z.edu + x.sex * y.att * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                           3.21888    0.20000  16.094  < 2e-16 ***
## x.sexLaki-laki                        0.18232    0.27080   0.673  0.50078    
## y.attYes                              1.60944    0.21909   7.346 2.04e-13 ***
## z.eduSarjana                          0.78846    0.24121   3.269  0.00108 ** 
## z.eduSMA                              1.18784    0.22846   5.199 2.00e-07 ***
## x.sexLaki-laki:y.attYes              -0.22314    0.29944  -0.745  0.45615    
## x.sexLaki-laki:z.eduSarjana          -0.09531    0.32891  -0.290  0.77199    
## x.sexLaki-laki:z.eduSMA              -0.27155    0.31442  -0.864  0.38778    
## y.attYes:z.eduSarjana                -0.71150    0.27127  -2.623  0.00872 ** 
## y.attYes:z.eduSMA                    -1.65945    0.27021  -6.141 8.19e-10 ***
## x.sexLaki-laki:y.attYes:z.eduSarjana  0.17250    0.37291   0.463  0.64367    
## x.sexLaki-laki:y.attYes:z.eduSMA      0.39832    0.37387   1.065  0.28670    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  2.1227e+02  on 11  degrees of freedom
## Residual deviance: -5.9952e-15  on  0  degrees of freedom
## AIC: 97.794
## 
## Number of Fisher Scoring iterations: 3

exp(model_saturated$coefficients)

##                          (Intercept)                       x.sexLaki-laki 
##                           25.0000000                            1.2000000 
##                             y.attYes                         z.eduSarjana 
##                            5.0000000                            2.2000000 
##                             z.eduSMA              x.sexLaki-laki:y.attYes 
##                            3.2800000                            0.8000000 
##          x.sexLaki-laki:z.eduSarjana              x.sexLaki-laki:z.eduSMA 
##                            0.9090909                            0.7621951 
##                y.attYes:z.eduSarjana                    y.attYes:z.eduSMA 
##                            0.4909091                            0.1902439 
## x.sexLaki-laki:y.attYes:z.eduSarjana     x.sexLaki-laki:y.attYes:z.eduSMA 
##                            1.1882716                            1.4893162

Interpretasi :

Koefisien, rata-rata log jumlah kasus untuk kategori referensi adalah 3,22
x.sexLaki-laki, Laki-laki tidak memiliki perbedaan yang signifikan dengan perempuan (non-signifikan, p-value = 0,5)
y.attYes, yang setuju akan penggunakan transportasi umum memiliki expected count sebesar 5 kali daripada yang tidak setuju (signifikan, p-value = 2.04e-16)
z.eduSarjana, Kelompok Sarjana memiliki expected count sebesar 2,2 kali dari kelompok Pascasarjana (signifikan, p-value = 0,001)
z.eduSMA, Kelompok SMA memiliki expected count sebesar 3,8 kali dari kelompok Pascasarjana (signifikan, p-value = 2e-07)
Interaksi Dua Arah dan Tiga Arah lainnya, sebagian besar tidak signifikan (p >0,05) kecuali y.attYes & z.eduSarjana

Model Homogenous

Model log-linear homogenous memasukkan semua efek utama dan semua interaksi dua arah, tanpa interaks 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} \]

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

## 
## Call:
## glm(formula = counts ~ x.sex + y.att + z.edu + x.sex * y.att + 
##     x.sex * z.edu + y.att * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  3.3139826  0.1573418  21.062  < 2e-16 ***
## x.sexLaki-laki               0.0004068  0.1621369   0.003 0.997998    
## y.attYes                     1.4941741  0.1646845   9.073  < 2e-16 ***
## z.eduSarjana                 0.7116467  0.1815336   3.920 8.85e-05 ***
## z.eduSMA                     1.0489892  0.1773990   5.913 3.36e-09 ***
## x.sexLaki-laki:y.attYes     -0.0004981  0.1393720  -0.004 0.997149    
## x.sexLaki-laki:z.eduSarjana  0.0512377  0.1543960   0.332 0.739995    
## x.sexLaki-laki:z.eduSMA     -0.0001531  0.1663356  -0.001 0.999266    
## y.attYes:z.eduSarjana       -0.6220797  0.1860067  -3.344 0.000825 ***
## y.attYes:z.eduSMA           -1.4564206  0.1864610  -7.811 5.68e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 212.2727  on 11  degrees of freedom
## Residual deviance:   1.2221  on  2  degrees of freedom
## AIC: 95.016
## 
## Number of Fisher Scoring iterations: 3

Apakah Terdapat Interaksi Tiga Arah (Saturated vs Homogenous)

Hipotesis

H0 : Tidak terdaat interaksi tiga arah (Model Homogenous sudah cukup)
H1 : Terdapat interaksi tiga arah(Model saturated diperlukan)

Hitung Selisih Deviace

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

## [1] 1.222056

Hitung Derajat Bebas

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

## [1] 2

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"

Interpretasi :

Pada taraf signifikansi 5%, belum ada cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi tiga arah antara jenis kelamin, sikap, dan tingkat pendidikan

Uji Model Interaksi Dua Arah (Homogenous vs Conditional on X)

Model Conditional on X

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

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

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

## 
## Call:
## glm(formula = counts ~ x.sex + y.att + z.edu + x.sex * y.att + 
##     x.sex * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  3.884e+00  1.041e-01  37.302  < 2e-16 ***
## x.sexLaki-laki              -1.453e-03  1.470e-01  -0.010   0.9921    
## y.attYes                     7.354e-01  9.556e-02   7.696  1.4e-14 ***
## z.eduSarjana                 2.364e-01  1.092e-01   2.164   0.0304 *  
## z.eduSMA                     6.454e-02  1.137e-01   0.568   0.5701    
## x.sexLaki-laki:y.attYes      2.149e-03  1.345e-01   0.016   0.9873    
## x.sexLaki-laki:z.eduSarjana  5.129e-02  1.536e-01   0.334   0.7384    
## x.sexLaki-laki:z.eduSMA      1.114e-11  1.607e-01   0.000   1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 212.273  on 11  degrees of freedom
## Residual deviance:  70.777  on  4  degrees of freedom
## AIC: 160.57
## 
## Number of Fisher Scoring iterations: 4

Pengujian Ada Tidaknya Interaksi Antara Y dan Z (Homogenous Model vs Conditional Association on X)

Hipotesis

H0 : Tidak ada interaksi antara Attitude (Y) dan Education (Z)
H1 : Ada interaksi antara Attitude (Y) dan Education (Z)

Taraf Signifikansi

α = 5%

Statistik Uji

Deviance.modelx <- model_conditional_X$deviance- model_homogenous$deviance
Deviance.modelx

## [1] 69.5552

Daerah Kritis

# Chi Square tabel dengan alpha = 0.05
 derajat.bebas <- (4- 2)
 derajat.bebas

## [1] 2

Chi Square Tabel

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

## [1] 5.991465

Keputusan

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

## [1] "Tolak"

Kesimpulan

Dengan taraf nyata 5%, terdapat cukup bukti untuk menolak 𝐻0 atau dapat dikatakan bahwa ada interaksi antara Attitude (Y) dan Education (Z). Model yang terbentuk tidak hanya sampai dua interaksi dengan X (conditional on X) sehingga interaksi Y*Z signifikan secara statistik

Uji Model Interaksi Dua Arah (Hoogenous vs Conditional on Y)

Model Conditional on Y

Model log-linear conditional pada Y memasukkan efek utama dan interaksi dua arah antara X dengan Y dan 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}_{jk} \]

model_conditional_Y <- glm(counts ~ x.sex + y.att + z.edu +
x.sex*y.att + y.att*z.edu,family = poisson(link = "log"))
summary(model_conditional_Y)

## 
## Call:
## glm(formula = counts ~ x.sex + y.att + z.edu + x.sex * y.att + 
##     y.att * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              3.304969   0.145933  22.647  < 2e-16 ***
## x.sexLaki-laki           0.018349   0.110605   0.166 0.868237    
## y.attYes                 1.492840   0.163928   9.107  < 2e-16 ***
## z.eduSarjana             0.737599   0.163943   4.499 6.82e-06 ***
## z.eduSMA                 1.048913   0.156688   6.694 2.17e-11 ***
## x.sexLaki-laki:y.attYes  0.002149   0.134501   0.016 0.987250    
## y.attYes:z.eduSarjana   -0.622086   0.185998  -3.345 0.000824 ***
## y.attYes:z.eduSMA       -1.456421   0.186461  -7.811 5.68e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 212.2727  on 11  degrees of freedom
## Residual deviance:   1.3792  on  4  degrees of freedom
## AIC: 91.173
## 
## Number of Fisher Scoring iterations: 3

Pengujian Ada Tidaknya Interaksi Antara X dan Z (Homogenous Model vs Conditional Association on Y)

Hipotesis

H0 : Tidak ada interaksi antara Jenis Kelamin (X) dan Education (Z)
H1 : Ada interaksi antara Jenis Kelamin (X) dan Education (Z)

Taraf Signifikansi

α = 5%

Statistik Uji

Deviance.modely <- model_conditional_Y$deviance- model_homogenous$deviance 
Deviance.modely

## [1] 0.1571842

Daerah Kritis

derajat.bebas <- (4- 2)
derajat.bebas

## [1] 2

Chi Square Tabel

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

## [1] 5.991465

Keputusan

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

## [1] "Terima"

Kesimpulan

Dengan taraf nyata 5%, tidak terdapat cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi antara Jenis Kelamin (X) dan Education (Z). Model yang terbentuk hanya sampai dua interaksi dengan Y (conditional on Y) sehingga interaksi X*Z tidak signifikan secara statistik.

Uji Model Interaksi Dua arah (Homogenous vs Conditional on Z)

Model Conditional on Z

Model log-linear conditional pada Z memasukkan efek utama dan interaksi dua arah antara X dengan 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^{XZ}_{ik} + \lambda^{YZ}_{jk} \]

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

## 
## Call:
## glm(formula = counts ~ x.sex + y.att + z.edu + x.sex * z.edu + 
##     y.att * z.edu, family = poisson(link = "log"))
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  3.314e+00  1.467e-01  22.595  < 2e-16 ***
## x.sexLaki-laki              -1.168e-13  1.155e-01   0.000 1.000000    
## y.attYes                     1.494e+00  1.492e-01  10.012  < 2e-16 ***
## z.eduSarjana                 7.116e-01  1.814e-01   3.923 8.76e-05 ***
## z.eduSMA                     1.049e+00  1.761e-01   5.957 2.58e-09 ***
## x.sexLaki-laki:z.eduSarjana  5.129e-02  1.536e-01   0.334 0.738443    
## x.sexLaki-laki:z.eduSMA      1.165e-13  1.607e-01   0.000 1.000000    
## y.attYes:z.eduSarjana       -6.221e-01  1.860e-01  -3.345 0.000824 ***
## y.attYes:z.eduSMA           -1.456e+00  1.865e-01  -7.811 5.68e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 212.2727  on 11  degrees of freedom
## Residual deviance:   1.2221  on  3  degrees of freedom
## AIC: 93.016
## 
## Number of Fisher Scoring iterations: 3

Pengujian Ada Tidaknya Interaksi Antara X dan Y (Homogenous Model vs Conditional Association on Z)

Hipotesis

H0 : Tidak ada interaksi antara Jenis Kelamin (X) dan Attitude (Y)
H1 : Ada interaksi antara Jenis Kelamin (X) dan Attitude (Y)

Taraf Signifikansi

α = 5%

Statistik Uji

Deviance.modelz <- model_conditional_Z$deviance- model_homogenous$deviance 
Deviance.modelz

## [1] 1.277092e-05

Daerah Kritis

derajat.bebas <- (3- 2)
derajat.bebas

## [1] 1

Chi Square Tabel

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

## [1] 3.841459

Keputusan

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

## [1] "Terima"

Kesimpulan

Dengan taraf nyata 5%, tidak terdapat cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi antara Jenis Kelamin (X) dan Attitude (Y). Model yang terbentuk hanya sampai dua interaksi dengan Z (conditional on Z) sehingga interaksi X*Y tidak signifikan secara statistik.

Model Terbaik

Dari seluruh pengujian, model terbaik adalah model yang signifikan yaitu model yang hanya memiliki interaksi antara Attitude (Y) dan Education (Z)

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

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

## 
## Call:
## glm(formula = counts ~ x.sex + y.att + z.edu + y.att * z.edu, 
##     family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            3.30424    0.13853  23.851  < 2e-16 ***
## x.sexLaki-laki         0.01980    0.06293   0.315 0.753025    
## y.attYes               1.49393    0.14921  10.012  < 2e-16 ***
## z.eduSarjana           0.73760    0.16394   4.499 6.82e-06 ***
## z.eduSMA               1.04891    0.15669   6.694 2.17e-11 ***
## y.attYes:z.eduSarjana -0.62209    0.18600  -3.345 0.000824 ***
## y.attYes:z.eduSMA     -1.45642    0.18646  -7.811 5.68e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 212.2727  on 11  degrees of freedom
## Residual deviance:   1.3795  on  5  degrees of freedom
## AIC: 89.174
## 
## Number of Fisher Scoring iterations: 3

Interpretasi Model Terbaik

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

##                              koef   exp_koef
## (Intercept)            3.30423567 27.2277228
## x.sexLaki-laki         0.01980263  1.0200000
## y.attYes               1.49392503  4.4545455
## z.eduSarjana           0.73759894  2.0909091
## z.eduSMA               1.04891262  2.8545455
## y.attYes:z.eduSarjana -0.62208606  0.5368234
## y.attYes:z.eduSMA     -1.45642063  0.2330690

Interpretasi :

x.sexLaki-laki, Tanpa memerhatikan sikap dan tingkat pendidikan, peluang seseorang berjenis kelamin laki-laki tidak berbeda signifikan dibanding perempuan sebesar 1,02 kali
y.attYes, Tanpa memerhatikan jenis kelamin dan tingkat pendidikan, peluang seseorang setuju akan penggunakan transportasi umum sebesar 4,45 kali daripada yang menolak penggunaan transportasi umum
z.eduSarjana, Tanpa memerhatikan jenis kelamin dan sikap, peluang seseorang memiliki tingkat pendidikan Sarjana yaitu sebesar 2,09 kali dibandingkan Pasca Sarjana
z.eduSMA, Tanpa memerhatikan jenis kelamin dan sikap, peluang seseorang memiliki tingkat pendidikan SMA yaitu sebesar 2,85 kali dibandingkan Pasca Sarjana
y.attYes & z.eduSarjana, Tanpa memerhatikan jenis kelamin, odds setuju akan penggunaan transportasi umum jika dia seorang Sarjana lebih rendah 46,4% daripada yang menolak
y.attYes & z.eduSMA, Tanpa memerhatikan jenis kelamin, odds setuju akan penggunaan transportasi umum jika dia seorang Sarjana lebih rendah 76,7% daripada yang menolak

Referensi

Fisher, R. A. (1935). The Design of Experiments. Oliver and Boyd.

Agresti, A. (2013). Categorical Data Analysis. John Wiley & Sons.

Mehta, C. R., & Patel, N. R. (1983). A Network Algorithm for Performing Fisher’s Exact Test in r × c Contingency Tables. Journal of the American Statistical Association, 78(382), 427-434.

Howell, D. C. (2012). Statistical Methods for Psychology (8th ed.). Cengage Learning

Anderson, C. J. (2018). Log‑linear Models for Contingency Tables

Fox, J. (2008). Applied Regression Analysis and Generalized Linear Models.

Christensen, R. (1997). Log-Linear Models and Logistic Regression. Springer.

Modul ADK (2025) Universitas Padjadjaran

Booklet Analisis Data Kategori

Maulana Hardy Rayyan

2025-04-13

Kata Pengantar

Pendahuluan

Tujuan Analisis Data Kategori

Manfaat Analisis Data Kategori

Metode dalam Analisis Data Kategori

Tabel Kontingensi dan Uji Chi-Square

Regresi Logistik

Analisis Korespondensi

Decision Tree dan Random Forest

Distribusi Probabilitas dalam Data Kategori

Distribusi Bernoulli

Distribusi Binomial

Distribusi Multinomial

Distribusi Poisson

Desain Sampling dalam Analisis Data Kategori

Eksperimen

Studi Cohort

Studi Cohort Prospective

Studi Cohort Retrospective

Studi Kasus Kontrol

Tabel Kontingensi 2 x 2

Data untuk Tabel Kontingensi 2 x 2

Distribusi Peluang dalam Tabel Kontingensi 2 x 2

Peluang Bersama

Peluang Marjinal

Peluang Bersyarat

Perhitungan menggunakan R

Ukuran Asosiasi dalam Dat Kategori 2x2

Risk Difference

Relative Risk

Odds Ratio

Perhitungan Menggunakan R

Inferensi Tabel Kontingensi Dua Arah

Estimasi

Estimasi Titik

Estimasi Interval

Perhitungan dengan R

Pengujian Hipotesis

Uji Proporsi

Proporsi Sampel

Proporsi Gabungan

Statistik Uji Z

Perhitungan menggunakan R

Uji Asosiasi

Risk Difference

Relative Risk

Odds ratio

Perhitungan Menggunakan R

Uji Independensi

Uji Chi-Square

Uji Likelihood Ratio

Uji Exact Fisher

Analisis Residual dalam Tabel Kontingensi

Jenis Residual

Tabel Kontingensi Tiga Arah

Tabel Parsial dan Marginal

Data untuk Tabel Kontingensi 3 arah

Freq Parsial Berdasarkan Lama Pendidikan

Distribusi Peluang

Peluang Bersama

Peluang Marginal

Peluang Bersyarat

Tabel Peluang Bersyarat

Ukuran Asosiasi

Risk Difference

Relative Risk

Odds Ratio

Inferensi

Uji CMH

Uji Homogenitas Breslow-Day

Generalized Linear Model

Exponential Family

Regresi Logistik

Regresi Poisson

Inferensi Generalized Linear Model

Ekspektasi dan Varians GLM

Metode Penaksiran Parameter