Analisis Data Kategori

Tabel Kontingensi 2x2

Tabel kontingensi 2×2 adalah tabel yang menyajikan distribusi frekuensi dua variabel kategorikal, masing-masing dengan dua kategori. Tabel ini membantu menganalisis hubungan atau asosiasi antara kedua variabel.

Notasi Tabel Kontingensi 2×2

Dalam analisis tabel kontingensi 2×2, kita menggunakan notasi berikut:

• \(n_{11}\) = Jumlah kasus positif Grup 1
  
• \(n_{12}\) = Jumlah kasus negatif Grup 1
  
• \(n_{21}\) = Jumlah kasus positif Grup 2
  
• \(n_{22}\) = Jumlah kasus negatif Grup 2

Total Baris:
- \(n_{1.} = n_{11} + n_{12}\) (Total Grup 1)
- \(n_{2.} = n_{21} + n_{22}\) (Total Grup 2)

Total Kolom:
- \(n_{.1} = n_{11} + n_{21}\) (Total Kategori Positif)
- \(n_{.2} = n_{12} + n_{22}\) (Total Kategori Negatif)

Grand Total:
- \(n = n_{11} + n_{12} + n_{21} + n_{22}\) (Total sampel)

Contoh Struktur Tabel:

	Positif (\(n_{.1}\))	Negatif (\(n_{.2}\))	Total
Grup 1	\(n_{11}\)	\(n_{12}\)	\(n_{1.}\)
Grup 2	\(n_{21}\)	\(n_{22}\)	\(n_{2.}\)
Total	\(n_{.1}\)	\(n_{.2}\)	\(n\)

Kasus :

Sebuah penelitian dilakukan untuk mengevaluasi hubungan antara kebiasaan minum kopi dan kejadian insomnia di kalangan pekerja kantoran. Peneliti mengumpulkan data dari 100 pekerja kantoran yang dibagi menjadi dua kelompok berdasarkan kebiasaan minum kopi mereka: kelompok yang minum kopi (K) dan kelompok yang tidak minum kopi (Kc). Selain itu, peneliti juga mencatat apakah pekerja tersebut mengalami insomnia (I) atau tidak mengalami insomnia (Ic). Berikut adalah hasil pengamatan dalam tabel kontingensi:

	Insomnia (I)	Tidak Insomnia (Ic)	Total
Minum Kopi (K)	30	20	50
Tidak Minum Kopi (Kc)	10	40	50
Total	40	60	100

Distribusi Peluang dalam tabel kontingensi 2x2

Peluang Bersama

Peluang bersama dihitung dengan rumus: \[ P(A \cap B) = \frac{n(A \cap B)}{n(S)} \]

# Menghitung peluang bersama
P_K_I <- 30/100  # P(K ∩ I)
P_K_Ic <- 20/100 # P(K ∩ Ic)
P_Kc_I <- 10/100 # P(Kc ∩ I)
P_Kc_Ic <- 40/100 # P(Kc ∩ Ic)

P_K_I; P_K_Ic; P_Kc_I; P_Kc_Ic

## [1] 0.3

## [1] 0.2

## [1] 0.1

## [1] 0.4

Interpretasi: Peluang bersama menunjukkan probabilitas terjadinya dua kejadian sekaligus, misalnya peluang seseorang yang minum kopi mengalami insomnia adalah 30%.

Peluang Marginal

Peluang marginal dihitung sebagai: \[ P(A) = \frac{n(A)}{n(S)} \]

# Menghitung peluang marginal
P_K <- (30 + 20)/100  # P(K)
P_Kc <- (10 + 40)/100 # P(Kc)
P_I <- (30 + 10)/100  # P(I)
P_Ic <- (20 + 40)/100 # P(Ic)

P_K; P_Kc; P_I; P_Ic

## [1] 0.5

## [1] 0.5

## [1] 0.4

## [1] 0.6

Interpretasi: Peluang marginal menunjukkan probabilitas terjadinya suatu kejadian tanpa mempertimbangkan kejadian lainnya, misalnya peluang seseorang minum kopi adalah 50%.

Peluang Bersyarat

Peluang bersyarat dihitung dengan rumus: \[ P(A | B) = \frac{P(A \cap B)}{P(B)} \]

# Menghitung peluang bersyarat
P_I_given_K <- P_K_I / P_K  # P(I | K)
P_Ic_given_K <- P_K_Ic / P_K # P(Ic | K)
P_I_given_Kc <- P_Kc_I / P_Kc # P(I | Kc)
P_Ic_given_Kc <- P_Kc_Ic / P_Kc # P(Ic | Kc)

P_I_given_K; P_Ic_given_K; P_I_given_Kc; P_Ic_given_Kc

## [1] 0.6

## [1] 0.4

## [1] 0.2

## [1] 0.8

Interpretasi: Peluang bersyarat menunjukkan probabilitas suatu kejadian terjadi jika diketahui bahwa kejadian lain sudah terjadi, misalnya peluang seseorang mengalami insomnia jika ia minum kopi adalah 60%.

Koefisien Phi,Cramer v, Koefisien lambda

Koefisien Phi (Metode Determinan)

\[ \phi = \frac{(ad - bc)}{\sqrt{(a+b)(c+d)(a+c)(b+d)}} \]

# Menghitung Koefisien Phi dengan metode determinan 
phi <- ((30 * 40) - (20 * 10)) / sqrt((30 + 20) * (10 + 40) * (30 + 10) * (20 + 40)) 
phi

## [1] 0.4082483

Interpretasi: Jika \(\phi\) mendekati 1 atau -1, maka hubungan antara variabel semakin kuat.

Koefisien Phi (Metode Chi-Square)

\[ \phi = \sqrt{\frac{X^2}{N}} \]

# Menghitung Koefisien Phi dengan metode Chi-Square 
X2 <- ((30*40 - 20*10)^2)*100/(50 * 50 * 40 * 60) 
n <- 100 # Total sampel 
phi_chi <- sqrt(X2 / n) 
phi_chi

## [1] 0.4082483

Interpretasi: Nilai \(\phi\) yang lebih tinggi menunjukkan hubunganyang lebih kuat.

Cramér’s V

\[V = \sqrt{\frac{X^2}{n(k-1)}}\]

# Menghitung Cramér's V 
k <- min(2,2) # Untuk tabel 2x2, k-1 = 1 
cramer_v <- sqrt(X2 / (n * (k - 1))) 
cramer_v

## [1] 0.4082483

Interpretasi: Cramér’s V memberikan ukuran kekuatan hubungan antara dua variabel kategorikal. Nilai mendekati 1 menunjukkan hubungan yang kuat.

Koefisien Lambda

\[ \lambda = \frac{E1 - E2}{E1} \]

# Menghitung Koefisien Lambda 
E1 <- max(30 + 20, 10 + 40) # Maksimal baris 
E2 <- 30 + 40 # Maksimum kolom
lambda <- (E1 - E2) / E1 
lambda

## [1] -0.4

Interpretasi: Nilai lambda menunjukkan seberapa besar pengurangan kesalahan dalam memprediksi insomnia berdasarkan konsumsi kopi. Jika mendekati 1, berarti variabel memiliki hubungan yang kuat.

Ukuran Asosiasi

Rasio Peluang (Odds Ratio, OR)

Odds Ratio membandingkan peluang terjadinya outcome antara kelompok terpapar dan tidak terpapar.

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

Keterangan:
- \(a\): Kasus pada terpapar
- \(b\): Bukan kasus pada terpapar
- \(c\): Kasus pada tidak terpapar
- \(d\): Bukan kasus pada tidak terpapar

Interpretasi:
- \(OR = 1\): Tidak ada asosiasi
- \(OR > 1\): Faktor risiko
- \(OR < 1\): Faktor protektif

Risiko Relatif (Relative Risk, RR)

Relative Risk membandingkan probabilitas outcome antara dua kelompok.

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

Interpretasi:
- \(RR = 1\): Risiko sama
- \(RR > 1\): Risiko meningkat
- \(RR < 1\): Risiko menurun

Selisih Risiko (Risk Difference, RD)

Mengukur perbedaan absolut risiko antara kelompok.

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

Interpretasi:
- \(RD = 0\): Tidak ada perbedaan
- \(RD > 0\): Risiko absolut lebih tinggi
- \(RD < 0\): Risiko absolut lebih rendah

Contoh pada kasus Insomnia

Odds Ratio :

\[ OR = \frac{(a/c)}{(b/d)} = \frac{(30/10)}{(20/40)} \]

# Menghitung Odds Ratio} 
OR <- (30/10) / (20/40) 
OR

## [1] 6

Interpretasi: Jika OR > 1, maka minum kopi meningkatkan risiko insomnia. Dalam kasus ini, OR > 1 menunjukkan adanya hubungan positif antara konsumsi kopi dan insomnia.

Risiko Relatif (Relative Risk) :

\[ RR = \frac{P(I | K)}{P(I | Kc)} \]

# Menghitung Risiko Relatif} 
RR <- P_I_given_K / P_I_given_Kc 
RR

## [1] 3

Interpretasi: Jika RR > 1, maka risiko insomnia lebih besar bagi peminum kopi dibandingkan yang tidak minum kopi.

Risk Difference (RD) :

n11 <- 30; n12 <- 20; n21<- 10 ; n22 <- 40
n1. <- n11 + n12; n2.<-n21+n22
p1<-(n11/n1.)
p2<-(n21/n2.)
rd <- p1 - p2
se_rd <- sqrt((p1 * (1 - p1) / n1.) + p2*((1 - p2) / n2.))
z_rd <- rd / se_rd
list(RD = rd, SE_RD = se_rd, Z_RD = z_rd)

## $RD
## [1] 0.4
## 
## $SE_RD
## [1] 0.08944272
## 
## $Z_RD
## [1] 4.472136

Inferensi Tabel Kontingensi 2 Arah

Estimasi

Estimasi bertujuan untuk memperkirakan parameter populasi berdasarkan data sampel. Estimasi dibagi menjadi:

Estimasi Titik

Digunakan untuk menentukan satu nilai spesifik sebagai perkiraan terbaik parameter populasi:

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

dimana:
- \(\hat{p}\): estimasi titik proporsi
- \(x\): jumlah “sukses” dalam sampel
- \(n\): total sampel

Estimasi Interval

Memberikan rentang nilai yang diyakini mengandung parameter populasi dengan tingkat kepercayaan tertentu:

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

dimana:
- \(Z_{\alpha/2}\): nilai distribusi normal standar (1.96 untuk CI 95%)

---

Contoh Kasus: Estimasi Proporsi Pengguna Smartphone X

Latar Belakang:
Survei terhadap 200 mahasiswa menemukan 75 orang menggunakan smartphone merek X.

1. Estimasi Titik

\[ \hat{p} = \frac{75}{200} = 0.375 \ \text{(37.5\%)} \]

p_hat <- 75/200
p_hat

## [1] 0.375

2. Estimasi Interval (95% CI)

\[ 0.375 \pm 1.96 \sqrt{\frac{0.375 \times 0.625}{200}} = 0.375 \pm 0.067 \]

Hasil:
\(\text{CI 95\%} = [0.308, 0.442]\) atau 30.8% hingga 44.2%

margin_error <- 1.96 * sqrt(p_hat*(1-p_hat)/200)
c(p_hat - margin_error, p_hat + margin_error)

## [1] 0.307904 0.442096

3. Visualisasi

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

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

Interpretasi

• Estimasi titik: 37.5% mahasiswa menggunakan smartphone X.
  
• Interval 95%: Proporsi sebenarnya di populasi diperkirakan antara 30.8% dan 44.2% dengan kepercayaan 95%.

Uji Hipotesis

1. Uji Proporsi Dua Sampel

Digunakan untuk membandingkan proporsi antara dua kelompok independen.

Hipotesis:

- \(H_0: p_1 = p_2\) (tidak ada perbedaan proporsi)

- \(H_1: p_1 \neq p_2\) (ada perbedaan proporsi)

Statistik Uji: \[ Z = \frac{\hat{p}_1 - \hat{p}_2}{\sqrt{\hat{p}(1-\hat{p})\left(\frac{1}{n_1} + \frac{1}{n_2}\right)}} \] di mana \(\hat{p} = \frac{x_1 + x_2}{n_1 + n_2}\) adalah proporsi gabungan.

Keputusan: Tolak \(H_0\) jika \(|Z| \geq z_{1-\alpha/2}\) atau p-value \(< \alpha\).

2. Uji Chi-Square Pearson

Digunakan untuk menguji independensi antara dua variabel kategorikal.

Hipotesis: - \(H_0\): Variabel independen - \(H_1\): Variabel dependen

Statistik Uji: \[ \chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \] dengan \(E_{ij} = \frac{(row\ total) \times (column\ total)}{grand\ total}\).

Asumsi:

- Sel yang diharapkan \(E_{ij} \geq 5\).

3. Uji Likelihood Ratio (G-Test)

Alternatif untuk uji chi-square ketika sampel kecil.

Statistik Uji: \[ G = 2 \sum O_{ij} \ln\left(\frac{O_{ij}}{E_{ij}}\right) \]

4. Uji Fisher Exact

Digunakan ketika sampel kecil atau ada sel dengan \(E_{ij} < 5\).

Hipotesis: Sama dengan chi-square.

Statistik Uji: \[ p = \frac{\binom{a+b}{a} \binom{c+d}{c}}{\binom{n}{a+c}} = \frac{(a+b)!(c+d)!(a+c)!(b+d)!}{a!b!c!d!n!} \]

5. Distribusi Hipergeometrik

Model probabilitas untuk tabel 2×2 dengan margin tetap.

Formula: \[ P(X = a) = \frac{\binom{a+b}{a} \binom{c+d}{c}}{\binom{n}{a+c}} \]

Catatan:

- Gunakan uji Fisher ketika sampel kecil (\(n < 20\))

- Uji chi-square dan G-test memerlukan \(E_{ij} \geq 5\)

- Distribusi hipergeometrik cocok untuk desain studi kasus-kontrol

Uji Proporsi

Proporsi Sampel

p1 <- 30 / 50 # proporsi Insomnia pada peminum kopi
p2 <- 10 / 50 # proporsi insomnia pada bukan peminum kopi 
p1; p2

## [1] 0.6

## [1] 0.2

Proporsi Gabungan

p1_2 <- (30 +10)/(50+50)
p1_2

## [1] 0.4

Statistik Uji Z

# Standard error
se <- sqrt(p1_2 * (1 - p1_2) * (1/50 + 1/50))

# Nilai Z
z <- (p1 - p2) / se

# Output
z

## [1] 4.082483

Interpretasi : Karena Z = 4.08 lebih besar dari nilai kritis untuk alpha = 0.05 (yaitu 1.96), kita menolak hipotesis nol, yang berarti ada perbedaan signifikan antara dua proporsi.

Uji Chi Square

set.seed(123)  
data <- matrix(c(30, 10, 20, 40), nrow = 2, byrow = TRUE)  
dimnames(data) <- list("Insomnia" = c("Ya", "Tidak"), "Kejadian" = c("Ya", "Tidak"))  
print(data)

##         Kejadian
## Insomnia Ya Tidak
##    Ya    30    10
##    Tidak 20    40

chisq_test <- chisq.test(data)  
print(chisq_test)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data
## X-squared = 15.042, df = 1, p-value = 0.0001052

Interpretasi : Dengan derajat bebas (df) = (2-1)(2-1) = 1, kita mendapatkan nilai Chi-Squared (15.042) dan karena P Value (0.0001052) < alpha (0,05) menyatakan bahwa ada hubungan signifikan antara variabel Insomnia dan Kopi.

Analisis Residual

Residual dalam tabel kontingensi digunakan untuk mengidentifkasi sel mana yang menyumbang paling banyak terhadap hubungan antara variabel kategori. Residual mengukur selisih antara frekuensi yang diamati dan frekuensi yang diharapkan berdasarkan model independensi. Jika suatu sel memiliki residual yang besar (positif atau negatif), berarti frekuensi observasi dalam sel tersebut sangat berbeda dari yang diharapkan. Sebaliknya, jika residualnya kecil, berarti nilai observasi mendekati nilai ekspektasi, sehingga sel tersebut tidak banyak menyumbang terhadap hubungan antara variabel Agresti, A. (2013).

# Data Observasi
observed <- matrix(c(30, 20, 10, 40), nrow = 2, byrow = TRUE)
# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected
# 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) * (1 - col_sum/total_sum) # Menampilkan hasil
list(
Pearson_Residual = pearson_residual,
Standardized_Residual = standardized_residual
)

## $Pearson_Residual
##           [,1]      [,2]
## [1,]  2.236068 -1.825742
## [2,] -2.236068  1.825742
## 
## $Standardized_Residual
##            [,1]       [,2]
## [1,]  0.6708204 -0.5477226
## [2,] -0.4472136  0.3651484

• Jika residual 0: – Tidak ada perbedaan signifkan antara jumlah observasi yang terjadi di suatu sel dengan jumlah yang diprediksi oleh model independensi. – Artinya, variabel baris dan kolom dalam tabel tidak memiliki hubungan kuat atau tidak menunjukkan pola keterkaitan yang jelas.

• Jika residual positif besar: – Frekuensi observasi jauh lebih tinggi dari yang diharapkan → menunjukkan hubungan positif yang kuat antara kategori tersebut.

• Jika residual negatif besar: – Frekuensi observasi jauh lebih rendah dari yang diharapkan → menunjukkan hubungan negatif atau tidak adanya asosiasi.

• Residual kecil (0) → Tidak ada hubungan antara variabel baris dan kolom.

• Residual besar (positif atau negatif) → Ada hubungan yang signifkan antara variabel baris dan kolom.

Mendeteksi Outlier menggunakan residual

# Data Observasi
observed <- matrix(c(30, 20, 10, 40), nrow = 2, byrow = TRUE)
colnames(observed) <- c("Insomnia", "tdk Insom")
rownames(observed) <- c("Ngopi", "")
# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected
# Pearson Residual
pearson_residual <- (observed - expected) / sqrt(expected)
42

## [1] 42

# 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) * (1 - col_sum/total_sum)# Menampilkan hasil
list(
Observed = observed,
Expected = expected,
Pearson_Residual = pearson_residual,
Standardized_Residual = standardized_residual
)

## $Observed
##       Insomnia tdk Insom
## Ngopi       30        20
##             10        40
## 
## $Expected
##       Insomnia tdk Insom
## Ngopi       20        30
##             20        30
## 
## $Pearson_Residual
##        Insomnia tdk Insom
## Ngopi  2.236068 -1.825742
##       -2.236068  1.825742
## 
## $Standardized_Residual
##         Insomnia  tdk Insom
## Ngopi  0.6708204 -0.5477226
##       -0.4472136  0.3651484

• Deteksi outlier menggunakan residual membantu mengidentifkasi kategori dalam tabel kontingensi yang sangat menyimpang dari ekspektasi.

• Jika |r| > 3, maka sel tersebut dapat dianggap sebagai outlier signifkan yang mungkin mempengaruhi kesimpulan analisis.

• Residual besar menunjukkan hubungan signifkan antara variabel, sedangkan residual mendekati nol menunjukkan tidak adanya hubungan.

Tabel kontingensi 3 arah

Tabel kontingensi 3 arah adalah representasi multidimensi dari distribusi frekuensi gabungan untuk tiga variabel kategorikal yang memungkinkan analisis interaksi kompleks. Tabel kontingensi tiga arah merupakan perpanjangan dari tabel kontingensi dua arah yang digunakan untuk menganalisis hubungan antara tiga variabel kategori secara simultan. Dalam banyak kasus, hubungan antara dua variabel (misalnya X dan Y) dapat dipengaruhi oleh variabel ketiga (Z), yang disebut sebagai variabel kontrol atau kovariat​

• Tabel Parsial: Tabel yang menyajikan hubungan antara X dan Y pada setiap kategori Z. Ini memungkinkan analisis hubungan bersyarat, di mana efek Z dikendalikan.

• Tabel Marginal: Tabel yang diperoleh dengan menjumlahkan semua kategori Z. Ini memberikan gambaran umum hubungan antara X dan Y tanpa mempertimbangkan Z, yang bisa menyebabkan distorsi interpretasi, seperti dalam Simpson’s Paradox

Dalam tabel kontingensi tiga arah :

• Risk Difference (RD) : ukuran absolut yang menghitung selisih probabilitas kejadian antara dua kelompok yang dibandingkan.

  Tujuan :

  • Menunjukkan dampak absolut dari faktor risiko terhadap probabilitas kejadian.

  • Memberikan informasi lebih konkret tentang berapa banyak kejadian yang dapat dicegah jika faktor risiko dihilangkan.

  • Berguna dalam pengambilan keputusan kebijakan kesehatan dan intervensi klinis.

  Interpretasi RD :

  • RD = 0 → Tidak ada perbedaan risiko antara kedua kelompok.

  • RD > 0 → Faktor risiko meningkatkan kemungkinan kejadian.

  • RD < 0 → Faktor risiko menurunkan kemungkinan kejadian.

  Rumus :

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

• Relative Risk (RR) atau risiko relatif : rasio probabilitas kejadian antara kelompok yang terpapar dengan kelompok yang tidak terpapar.

  • Menentukan seberapa besar risiko suatu kejadian pada kelompok yang terpapar dibandingkan dengan yang tidak terpapar.
  • Digunakan dalam studi kohort dan uji klinis.
  • Memberikan pemahaman tentang potensi faktor risiko terhadap suatu penyakit.

  Interpretasi RR :

  • RR = 1 → Tidak ada hubungan antara faktor risiko dan kejadian.
  • RR > 1 → Faktor risiko meningkatkan kemungkinan kejadian.
  • RR < 1 → Faktor risiko bersifat protektif (mengurangi kemungkinan kejadian).

  Rumus :

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

• Odds Ratio (OR) : perbandingan peluang kejadian antara dua kelompok, biasanya digunakan dalam studi kasus-kontrol.

  Tujuan :

  • Memberikan estimasi hubungan antara faktor risiko dan kejadian.
  • Digunakan dalam kondisi ketika tidak memungkinkan menghitung risiko absolut.
  • Berguna dalam epidemiologi dan penelitian klinis.

  Interpretasi OR :

  • OR = 1 → Tidak ada hubungan antara faktor risiko dan kejadian.
  • OR > 1 → Faktor risiko meningkatkan peluang kejadian.
  • OR < 1 → Faktor risiko mengurangi peluang kejadian.

  Rumus :

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

Konsep Matematis

Untuk tiga variabel \(X\), \(Y\), dan \(Z\) dengan level \(I\), \(J\), dan \(K\):

\[ \pi_{ijk} = P(X = i, Y = j, Z = k) \]

dengan batasan: \[ \sum_{i=1}^I \sum_{j=1}^J \sum_{k=1}^K \pi_{ijk} = 1 \]

Model Log-linear yang sering digunakan: \[ \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} \]

Tujuan Analisis

Menurut Fienberg (2007):

1. Menguji Interaksi Tiga Faktor:
   Contoh: Apakah efek merokok \((X)\) pada kanker \((Y)\) berbeda berdasarkan jenis kelamin \((Z)\)?

2. Analisis Asosiasi Bersyarat:
   \[ \theta_{XY(k)} = \frac{\pi_{11k} \pi_{22k}}{\pi_{12k} \pi_{21k}} \] Mengukur asosiasi \(X\) dan \(Y\) pada level tertentu \(Z=k\).
3. Identifikasi Variabel Pengecoh (Confounding):
   Memisahkan efek langsung dan tidak langsung antar variabel.

Kelebihan dan Keterbatasan

Kelebihan (Pearson, 1904):

- Dapat menangkap hubungan non-linear yang tidak terlihat di analisis 2 arah.

Keterbatasan (Cochran, 1954):

- Masalah Sparsitas: Banyak sel kosong jika sampel kecil.

- Interpretasi Kompleks: Interaksi orde tinggi sulit dijelaskan secara substantif.

Contoh Kasus

Data ini diperoleh dari survei kesehatan masyarakat di Indonesia, yang mencatat jumlah individu yang terdiagnosis AIDS (Yes) atau tidak terdiagnosis AIDS (No) berdasarkan gender dan orientasi seksual.

Gender	Orientasi Seksual	Yes	No
Laki-laki	Heteroseksual	50	100
	Homoseksual	80	70
Perempuan	Heteroseksual	90	30
	Homoseksual	20	60

Variabel yang digunakan

X : Orientasi Seksual (Heteroseksual / Homoseksual / Biseksual)

Y : Diagnosis Aids (Ya untuk Terdiagnosis AIDS, Tidak untuk Tidak Terdiagnosis AIDS)

Z : Gender (Laki laki / Perempuan)

Tabel Parsial Z₁

Orientasi Seksual	\(P(Y=\text{Yes} \mid X, Z=Laki-laki)\)	\(P(Y=\text{No} \mid X, Z=Laki-laki)\)
Heteroseksual	\(0.333\)	\(0.667\)
Homoseksual	\(0.533\)	\(0.467\)

Menghitung Risk Difference (RD) atau Beda Peluang

\[ RD = P(Y = \text{Yes} \mid X = \text{Homoseksual}, Z = \text{Laki-laki}) - P(Y = \text{Yes} \mid X = \text{Heteroseksual}, Z = \text{Laki-laki}) \]

\[ P(Y = \text{Yes} \mid X = \text{Heteroseksual}, Z = \text{Laki-laki}) = \frac{50}{150} = 0.333 \]

\[ P(Y = \text{Yes} \mid X = \text{Homoseksual}, Z = \text{Laki-laki}) = \frac{80}{150} = 0.533 \]

# Data 
hetero_yes <- 50 
hetero_total <- 150
homo_yes <- 80
homo_total <- 150
# Perhitungan peluang bersyarat
P_hetero <- hetero_yes / hetero_total
P_homo <- homo_yes / homo_total  
# Risk Difference (RD)
RD <- P_hetero - P_homo 
RD 

## [1] -0.2

\[ RD = 0.333 - 0.533 = -0.200 \]

Interpretasi :
Artinya gender laki laki yang homoseksual memiliki peluang lebih besar 2% untuk terkena AIDS daripada laki laki yang heteroseksual.

Menghitung Resiko Relatif (RR)

\[ RR = \frac{P(Y=Yes | X=1, Z=1)}{P(Y=Yes | X=0, Z=1)} \]

# Data jumlah kasus
yes_homo <- 80   # Yes pada laki-laki homoseksual
total_homo <- 150 # Total laki-laki homoseksual (Yes + No)
yes_hetero <- 50  # Yes pada laki-laki heteroseksual
total_hetero <- 150 # Total laki-laki heteroseksual (Yes + No)
# Menghitung peluang bersyarat
p_yes_homo <- yes_homo / total_homo
p_yes_hetero <- yes_hetero / total_hetero  
# Menghitung Risiko Relatif (RR) 
RR <- p_yes_homo / p_yes_hetero 
RR

## [1] 1.6

Dengan data yang diberikan:

\[ RR = \frac{0.533}{0.333} = 1.60 \]

Interpretasi :

Artinya Gender laki laki yang homoseksual yang terkena HIV atau AIDS 1,6 kali lebih banyak daripada laki laki yang hetero seksual.

Perhitungan Odds Ratio (OR) untuk Laki-laki

Rumus Odds Ratio (OR)

Odds Ratio dihitung menggunakan rumus:

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

Dimana:

a = Kasus (Yes) pada kelompok homoseksual = 80

b = Non-kasus (No) pada kelompok homoseksual = 70

c = Kasus (Yes) pada kelompok heteroseksual = 50

d = Non-kasus (No) pada kelompok heteroseksual = 100

Implementasi dalam R

# Data untuk laki-laki
a <- 80  # Yes pada homoseksual
b <- 70  # No pada homoseksual 
c <- 50  # Yes pada heteroseksual 
d <- 100 # No pada heteroseksual  
# Menghitung Odds Ratio
OR <- (a * d) / (b * c) 
OR

## [1] 2.285714

Interpretasi :

Berdasarkan perhitungan, OR = 2.29, yang berarti laki-laki homoseksual memiliki 2.29 kali lebih besar peluang untuk terdiagnosis AIDS dibanding laki-laki heteroseksual.

Tabel Parsial Z₂

Peluang bersyarat dihitung dengan rumus:

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

Di mana:

- \(Y\) adalah diagnosis AIDS (Yes/No).

- \(X\) adalah orientasi seksual (Heteroseksual/Homoseksual).

- \(Z\) adalah gender (Perempuan).

Peluang Bersyarat nya disajikan dalam tabel sebagai berikut :

Gender	Orientasi Seksual	P(Yes)	P(No)
Perempuan	Heteroseksual	0.750	0.250
Perempuan	Homoseksual	0.250	0.750

Menghitung Beda Peluang atau Risk Difference

\[ RD = P(Y = \text{Yes} \mid X = \text{Heteroseksual}, Z = \text{Perempuan}) - P(Y = \text{Yes} \mid X = \text{Homoseksual}, Z = \text{Perempuan}) \]

\[ P(Y = \text{Yes} \mid X = \text{Heteroseksual}, Z = \text{Perempuan}) = \frac{90}{120} = 0.750 \]

\[ P(Y = \text{Yes} \mid X = \text{Homoseksual}, Z = \text{Perempuan}) = \frac{20}{80} = 0.250 \]

# Data perempuan 
hetero_yes <- 90
hetero_total <- 120
homo_yes <- 20
homo_total <- 80  
# Perhitungan peluang bersyarat
P_hetero <- hetero_yes / hetero_total
P_homo <- homo_yes / homo_total
#Risk Difference (RD)
RD <- P_hetero - P_homo
RD 

## [1] 0.5

Interpretasi :

Artinya gender perempuan yang heteroseksual memiliki peluang lebih besar 12,5% untuk terkena HIV atau aids daripada perempuan yang homoseksual

Menghitung Resiko relatif

\[ P(Y = \text{Yes} \mid X = \text{Heteroseksual}, Z = \text{Perempuan}) = \frac{90}{120} = 0.750 \]

\[ P(Y = \text{Yes} \mid X = \text{Homoseksual}, Z = \text{Perempuan}) = \frac{20}{80} = 0.250 \]

# Data perempuan 
hetero_yes <- 90
hetero_total <- 120
homo_yes <- 20
homo_total <- 80
# Perhitungan peluang bersyarat
P_hetero <- hetero_yes / hetero_total
P_homo <- homo_yes / homo_total
# Relative Risk (RR)
RR <- P_homo / P_hetero 
RR 

## [1] 0.3333333

\[ RR = \frac{0.250}{0.750} = 0.333 \]

Interpretasi :

• RR < 1 menunjukkan bahwa perempuan homoseksual memiliki risiko lebih rendah dibandingkan perempuan heteroseksual.
• Nilai RR = 0.333 berarti perempuan homoseksual memiliki 33.3% kemungkinan terkena kondisi tersebut dibandingkan perempuan heteroseksual.

Menghitung Odds Ratio

\[ P(Y = \text{Yes} | X = \text{Heteroseksual}, Z = \text{Perempuan}) = \frac{90}{120} = 0.75 \]

\[ P(Y = \text{Yes} | X = \text{Homoseksual}, Z = \text{Perempuan}) = \frac{20}{80} = 0.250 \]

\[ \text{Odds Heteroseksual} = \frac{0.75}{1 - 0.75} = \frac{0.75}{0.25} = 3 \]

\[ \text{Odds Homoseksual} = \frac{0.250}{1 - 0.250} = \frac{0.250}{0.750} = 0,333 \]

\[ OR = \frac{3}{0.333} = 9 \]

# Data 
hetero_yes <- 90
hetero_total <- 120 # Jumlah total perempuan heteroseksual
homo_yes <- 20
homo_total <- 80  # Jumlah total perempuan homoseksual
# Perhitungan peluang bersyarat
P_hetero <- hetero_yes / hetero_total
P_homo <- homo_yes / homo_total
# Perhitungan Odds
odds_hetero <- P_hetero / (1 - P_hetero)
odds_homo <- P_homo / (1 - P_homo)
# Perhitungan Odds Ratio (OR) 
OR <- odds_hetero / odds_homo 
OR 

## [1] 9

Interpretasi :

Odds Ratio (OR) = 9, yang berarti bahwa perempuan heteroseksual memiliki 9 kali kemungkinan lebih besar mengalami kejadian Y = Yes dibandingkan perempuan homoseksual.

Tabel Marginal

Orientasi Seksual	Yes	No
Heteroseksual	140	130
Homoseksual	100	130

Peluang bersyaratnya disajikan dalam tabel berikut :

orientasi Seksual	P(Yes)	P(No)
Heteroseksual	0.518	0.482
Homoseksual	0.434	0.565

Menghitung Beda Peluang atau Risk Difference

\[ RD = P(Y = Yes | X = \text{Heteroseksual}) - P(Y = Yes | X = \text{Homoseksual}) \]

\[ RD = 0.518 - 0.434 = 0.08 \]

# Data
hetero_yes <- 140
hetero_total <- 140 + 130
homo_yes <- 100
homo_total <- 100 + 130 
# Perhitungan peluang bersyarat
P_hetero <- hetero_yes / hetero_total
P_homo <- homo_yes / homo_total  
# Risk Difference (RD)
RD <- P_hetero - P_homo 
RD 

## [1] 0.08373591

Interpretasi :
Artinya homoseksual memiliki peluang lebih kecil 8% untuk terkena AIDS daripada yang heteroseksual.

Menghitung Resiko Relatif

\[ RR = \frac{P(Y = Yes | X = \text{Heteroseksual})}{P(Y = Yes | X = \text{Homoseksual})} \]

\[ RR = \frac{0.518}{0.434} = 1.19 \]

# Data
hetero_yes <- 140
hetero_total <- 140 + 130
homo_yes <- 100
homo_total <- 100 + 130
# Perhitungan peluang bersyarat
P_hetero <- hetero_yes / hetero_total
P_homo <- homo_yes / homo_total 
# Relative Risk (RR) 
RR <- P_hetero / P_homo
RR 

## [1] 1.192593

Intepretasi :

• RR = 1.19 berarti kelompok hetero memiliki peluang 19% lebih tinggi untuk mengalami outcome (yes) dibanding kelompok homo.

• Contoh: Jika outcome adalah “sembuh dari penyakit”, kelompok hetero 1.19 kali lebih mungkin sembuh daripada kelompok homo.

• Ada asosiasi positif antara kelompok hetero dan outcome. Namun:

  RR mendekati 1 (1.19) menunjukkan efek yang lemah.

  Perlukan uji statistik (seperti p-value atau interval kepercayaan) untuk menentukan signifikansi.

  Catatan:

• Jika RR > 1: Kelompok pertama lebih berisiko/berpeluang.

• Jika RR = 1: Tidak ada perbedaan.

• Jika RR < 1: Kelompok pertama lebih rendah risikonya

Menghitung Odds Ratio (OR)

\[ OR = \frac{\frac{P(Y = \text{Yes} | X = \text{Heteroseksual})}{P(Y = \text{No} | X = \text{Heteroseksual})}}{\frac{P(Y = \text{Yes} | X = \text{Homoseksual})}{P(Y = \text{No} | X = \text{Homoseksual})}} \]

\[ OR = \frac{\frac{140}{130}}{\frac{100}{130}} \]

\[ OR = \frac{1.077}{0.769} = 1.4 \]

# Data jumlah kasus
hetero_yes <- 140
hetero_no <- 130
homo_yes <- 100 
homo_no <- 130  
# Menghitung odds untuk masing-masing kelompok
odds_hetero <- hetero_yes / hetero_no 
odds_homo <- homo_yes / homo_no
# Menghitung Odds Ratio (OR)
OR <- odds_hetero / odds_homo
OR 

## [1] 1.4

Interpretasi :

• Kelompok hetero memiliki odds (peluang relatif) 1.4 kali lebih tinggi untuk mengalami outcome (yes) dibanding kelompok homo.

• Dengan kata lain, kelompok hetero 40% lebih mungkin mengalami outcome daripada kelompok homo.

Conditional Dependence

Conditional Dependence adalah konsep dalam teori probabilitas dan statistika yang menggambarkan hubungan antara dua atau lebih variabel acak yang menjadi dependen (saling bergantung) ketika nilai dari satu atau lebih variabel lain diketahui. Conditional independence (kemandirian bersyarat) terjadi ketika dua variabel \(X\) dan \(Y\) menjadi independen setelah dikontrol oleh variabel ketiga \(Z\).

Tujuan

1. Mengungkap Hubungan Tersembunyi – Mengetahui ketergantungan antar variabel ketika faktor lain (ZZ) diperhitungkan.

2. Meningkatkan Akurasi Model – Memperbaiki prediksi dengan memanfaatkan struktur ketergantungan bersyarat.

3. Menghindari Bias Analisis – Mencegah kesimpulan salah dengan mengontrol variabel perancu (ZZ).

Definisi Formal

Variabel \(X\) dan \(Y\) bersyarat independen terhadap \(Z\) jika:

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

Representasi Frekuensi

Dalam notasi frekuensi tabel kontingensi:

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

dengan notasi:

- \(n_{ijk}\): Frekuensi sel untuk \(X=i\), \(Y=j\), \(Z=k\)
- \(n_{i+k}\): Total baris untuk \(X=i\) pada \(Z=k\)

- \(n_{+jk}\): Total kolom untuk \(Y=j\) pada \(Z=k\)

- \(n_{k++}\): Total marginal untuk \(Z=k\)

Contoh :

# Membuat tabel untuk laki-laki
data_laki <- matrix(c(50, 100, 80, 70), nrow = 2, byrow = TRUE)
rownames(data_laki) <- c("Heteroseksual", "Homoseksual")
colnames(data_laki) <- c("Yes", "No")

# Membuat tabel untuk perempuan
data_perempuan <- matrix(c(90, 30, 20, 60), nrow = 2, byrow = TRUE)
rownames(data_perempuan) <- c("Heteroseksual", "Homoseksual")
colnames(data_perempuan) <- c("Yes", "No")

# Menampilkan tabel
cat("Tabel untuk Laki-laki:\n")

## Tabel untuk Laki-laki:

print(data_laki)

##               Yes  No
## Heteroseksual  50 100
## Homoseksual    80  70

cat("\nTabel untuk Perempuan:\n")

## 
## Tabel untuk Perempuan:

print(data_perempuan)

##               Yes No
## Heteroseksual  90 30
## Homoseksual    20 60

# Uji Chi-square untuk conditional independence
cat("\nUji Chi-square untuk Laki-laki:\n")

## 
## Uji Chi-square untuk Laki-laki:

print(chisq.test(data_laki))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data_laki
## X-squared = 11.416, df = 1, p-value = 0.000728

cat("\nUji Chi-square untuk Perempuan:\n")

## 
## Uji Chi-square untuk Perempuan:

print(chisq.test(data_perempuan))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data_perempuan
## X-squared = 46.486, df = 1, p-value = 9.229e-12

# Uji Fisher Exact untuk perempuan karena ada sel dengan nilai <5
cat("\nUji Fisher Exact untuk Perempuan:\n")

## 
## Uji Fisher Exact untuk Perempuan:

print(fisher.test(data_perempuan))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  data_perempuan
## p-value = 2.844e-12
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   4.471007 18.294079
## sample estimates:
## odds ratio 
##   8.880715

# Uji Fisher Exact untuk laki laki  karena ada sel dengan nilai <5
cat("\nUji Fisher Exact untuk laki:\n")

## 
## Uji Fisher Exact untuk laki:

print(fisher.test(data_laki))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  data_laki
## p-value = 0.0007004
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.2667511 0.7164144
## sample estimates:
## odds ratio 
##  0.4387393

Interpretasi Hasil Uji Statistik

1. Uji Chi-square (\(\chi^2\))

• Untuk Laki-laki:
  • \(\chi^2 = 11.416\), \(p\text{-value} = 0.000728\)
  • Kesimpulan: Terdapat hubungan signifikan (\(p < 0.05\)) antara variabel.
  • Efek: Kuat (\(\chi^2 > 3.84\) untuk \(df = 1\)).
• Untuk Perempuan:
  • \(\chi^2 = 46.486\), \(p\text{-value} = 9.229 \times 10^{-12}\)
  • Kesimpulan: Hubungan sangat signifikan (\(p \ll 0.001\)).
  • Efek: Sangat kuat.

Formula Chi-square: \[ \chi^2 = \sum \frac{(O_i - E_i)^2}{E_i} \]

2. Uji Fisher Exact

• Untuk Perempuan:
  • \(OR = 8.881\), 95% CI [4.471, 18.294]
  • \(p\text{-value} = 2.844 \times 10^{-12}\)
  • Kesimpulan: Odds kelompok eksperimen 8.9x lebih tinggi (signifikan).
• Untuk Laki-laki:
  • \(OR = 0.439\), 95% CI [0.267, 0.716]
  • \(p\text{-value} = 0.0007004\)
  • Kesimpulan: Odds kelompok eksperimen 56% lebih rendah (signifikan).

Catatan:

- \(OR > 1\): Asosiasi positif

- \(OR < 1\): Asosiasi negatif

- CI tidak mencakup 1 \(\Rightarrow\) signifikan

Inferensi Tabel Kontingensi 3 Arah

Independensi Bersyarat dalam Tabel Kontingensi Tiga Arah

1. Konsep Independensi Bersyarat

Dua variabel \(X\) dan \(Y\) dikatakan independent bersyarat terhadap variabel ketiga \(Z\) jika: \[ OR(X,Y|Z) = 1 \] Artinya, setelah mengontrol \(Z\), tidak ada hubungan antara \(X\) dan \(Y\) dalam setiap strata \(Z\).

Catatan: - Fenomena \(X \perp Y\) di setiap strata \(Z\) tetapi tidak marginal disebut Paradoks Simpson. - Asumsi ini banyak digunakan dalam epidemiologi dan model machine learning.

---

2. Uji Cochran-Mantel-Haenszel (CMH)

Tujuan:
Menguji hubungan \(X\) dan \(Y\) dengan mengontrol confounder \(Z\) dalam tabel kontingensi berlapis.

Hipotesis: - \(H_0\): \(\theta_{xy(k)} = 1\) untuk semua strata \(k\)
(Tidak ada hubungan di setiap lapisan \(Z\)) - \(H_1\): \(\theta_{xy(k)} \neq 1\) untuk minimal satu strata \(k\)

Statistik Uji CMH: \[ CMH = \frac{\left[\sum_k (n_{11k} - \mu_{11k})\right]^2}{\sum_k \text{var}(n_{11k})} \] dengan: - \(\mu_{11k} = E(n_{11k}) = \frac{n_{1,k} \cdot n_{.1k}}{n_{..k}}\)
- \(\text{var}(n_{11k}) = \frac{n_{1,k} \cdot n_{2,k} \cdot n_{.1k} \cdot n_{.2k}}{n_{..k}^2 (n_{..k}-1)}\)

Keputusan: - Tolak \(H_0\) jika \(CMH > \chi^2_{1, \alpha}\) atau \(p\text{-value} < \alpha\).

---

# Data dalam bentuk array 3D (2 strata x 2x2)
data <- array(c(50, 80, 100, 70,   # Strata 1: Laki-laki
                90, 20, 30, 60),    # Strata 2: Perempuan
              dim = c(2, 2, 2),
              dimnames = list(
                Orientasi = c("Heteroseksual", "Homoseksual"),
                Outcome = c("Yes", "No"),
                Gender = c("Laki-laki", "Perempuan")
              ))

# Uji CMH
cmh_test <- mantelhaen.test(data)
print(cmh_test)

## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  data
## Mantel-Haenszel X-squared = 2.3752, df = 1, p-value = 0.1233
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.9309767 1.8247239
## sample estimates:
## common odds ratio 
##          1.303371

# Hitung OR per strata
or_laki <- (50 * 70) / (80 * 100)  # OR Laki-laki
or_perempuan <- (90 * 60) / (20 * 30)  # OR Perempuan
cat("\nOdds Ratio per Strata:",
    "\n- Laki-laki:", round(or_laki, 2),
    "\n- Perempuan:", round(or_perempuan, 2))

## 
## Odds Ratio per Strata: 
## - Laki-laki: 0.44 
## - Perempuan: 9

# Hitung OR Mantel-Haenszel manual
or_mh <- (50 * 70 / 300 + 90 * 60 / 200) / (80 * 100 / 300 + 20 * 30 / 200)
cat("\n\nOdds Ratio Gabungan (Mantel-Haenszel):", round(or_mh, 2))

## 
## 
## Odds Ratio Gabungan (Mantel-Haenszel): 1.3

Interpretasi Hasil Uji Cochran-Mantel-Haenszel (CMH)

1. Hasil Uji CMH

• Statistik Uji:
  \(\text{Mantel-Haenszel } \chi^2 = 2.3752\), \(df = 1\), \(p\text{-value} = 0.1233\)
  
• Kesimpulan:
  Karena \(p\text{-value} > 0.05\), gagal tolak \(H_0\). Tidak ada hubungan signifikan antara orientasi seksual dan outcome setelah mengontrol gender (\(\alpha = 0.05\)).
  
• Interval Kepercayaan OR Gabungan:
  95% CI = [0.931, 1.825] (mencakup 1 \(\Rightarrow\) tidak signifikan).

2. Odds Ratio (OR) per Strata

• Laki-laki:
  \(OR = 0.44\)
  Interpretasi:
  • Kelompok homoseksual memiliki odds 56% lebih rendah untuk outcome “Yes” dibanding heteroseksual pada strata laki-laki.
    
  • Catatan: \(\frac{1}{0.44} \approx 2.27\) → Heteroseksual 2.27x lebih mungkin “Yes”.
• Perempuan:
  \(OR = 9.0\)
  Interpretasi:
  • Kelompok homoseksual memiliki odds 9 kali lebih tinggi untuk outcome “Yes” dibanding heteroseksual pada strata perempuan.

3. OR Gabungan (Mantel-Haenszel)

• Estimasi: \(OR_{MH} = 1.30\)
  Interpretasi:
  • Secara gabungan (setelah dikontrol gender), tidak ada efek yang signifikan (karena CI mencakup 1).
    
  • Ketidakseragaman OR: Perbedaan besar antara strata (\(0.44\) vs \(9.0\)) menunjukkan interaksi antara gender dan orientasi seksual.

Uji Homogenitas Odds Ratio dengan Breslow-Day

1. Definisi Homogenitas

• Asosiasi homogen:
  \(\theta_{xy(1)} = \theta_{xy(2)} = \dots = \theta_{xy(K)}\)
  (OR sama di semua strata \(Z\)).
• Jika OR berbeda antar strata \(\Rightarrow\) interaksi antara \(X\), \(Y\), dan \(Z\).

2. Hipotesis

• \(H_0\): \(\theta_{XY(1)} = \theta_{XY(2)} = \dots = \theta_{XY(K)}\)
  
• \(H_1\): Minimal satu OR berbeda.

3. Statistik Uji Breslow-Day

\[ X_{\text{HBD}}^2 = \sum_{j=1}^K \frac{(a_j - \tilde{a}_j)^2}{\text{Var}(a_j|\hat{OR}_{MH})} \] Keterangan: - \(a_j\): Kasus terpapar diamati di strata \(j\). - \(\tilde{a}_j\): Kasus terpapar diharapkan di bawah \(H_0\). - \(\text{Var}(a_j|\hat{OR}_{MH})\): Varians kondisional.

4. Koreksi Tarone

\[ X_{\text{HBDT}}^2 = X_{\text{HBD}}^2 - \frac{\left(\sum_{j=1}^K (a_j - \tilde{a}_j)\right)^2}{\sum_{j=1}^K \text{Var}(a_j|\hat{OR}_{MH})} \]

5. Keputusan Hipotesis

• Derajat kebebasan: \(df = K - 1\)
  
• Tolak \(H_0\) jika:
  \(X_{\text{HBDT}}^2 > \chi^2_{\alpha, df}\) atau \(p\text{-value} < 0.05\).

6. Langkah Perhitungan

1. Hitung \(\hat{OR}_{MH}\): \[ \hat{OR}_{MH} = \frac{\sum_{j=1}^K \frac{a_j d_j}{n_j}}{\sum_{j=1}^K \frac{b_j c_j}{n_j}} \]
2. Cari \(\tilde{a}_j\) dari persamaan kuadrat**: \[ (1-\hat{OR}_{MH})\tilde{a}_j^2 + \left[n_{2j} - m_{1j} + \hat{OR}_{MH}(n_{1j} + m_{1j})\right]\tilde{a}_j - m_{1j}n_{1j}\hat{OR}_{MH} = 0 \]
3. Hitung varians: \[ \text{Var}(a_j|\hat{OR}_{MH}) = \left(\frac{1}{\tilde{a}_j} + \frac{1}{\tilde{b}_j} + \frac{1}{\tilde{c}_j} + \frac{1}{\tilde{d}_j}\right)^{-1} \]
4. Statistik uji + koreksi Tarone (Langkah 3-4).

7. Interpretasi

• \(p\text{-value} < 0.05\): Tolak \(H_0\) \(\Rightarrow\) OR tidak homogen (ada interaksi).
• \(p\text{-value} \geq 0.05\): Gagal tolak \(H_0\) \(\Rightarrow\) OR homogen.

Contoh Kasus Uji Homogenitas Breslow-Day

Studi Kasus: Pengaruh Merokok (\(X\)) terhadap Kanker Paru (\(Y\)) dengan Stratifikasi Usia (\(Z\)).

Data Tabel Kontingensi Berstrata

(Frekuensi kasus kanker paru pada perokok vs non-perokok, dibagi per kelompok usia)

Usia (Strata)	Merokok	Kanker: Ya	Kanker: Tidak	Total
<40 tahun (\(j=1\))	Ya	20 (\(a_1\))	80 (\(b_1\))	100 (\(m_{1}\))
	Tidak	10 (\(c_1\))	90 (\(d_1\))	100 (\(m_{2}\))
Total		30 (\(n_{1}\))	170 (\(n_{2}\))	200 (\(N_1\))
≥40 tahun (\(j=2\))	Ya	60 (\(a_2\))	40 (\(b_2\))	100 (\(m_{1}\))
	Tidak	30 (\(c_2\))	70 (\(d_2\))	100 (\(m_{2}\))
Total		90 (\(n_{1}\))	110 (\(n_{2}\))	200 (\(N_2\))

---

Langkah 1: Hitung Odds Ratio (OR) per Strata

• Strata <40 tahun:
  \[ OR_1 = \frac{a_1 \times d_1}{b_1 \times c_1} = \frac{20 \times 90}{80 \times 10} = 2.25 \]
  
• Strata ≥40 tahun:
  \[ OR_2 = \frac{60 \times 70}{40 \times 30} = 3.5 \]

Catatan: OR berbeda antar strata \(\Rightarrow\) Diduga ada interaksi usia.

---

Langkah 2: Estimasi OR Gabungan (Mantel-Haenszel)

\[ \hat{OR}_{MH} = \frac{\frac{20 \times 90}{200} + \frac{60 \times 70}{200}}{\frac{80 \times 10}{200} + \frac{40 \times 30}{200}} = \frac{9 + 21}{4 + 6} = 3.0 \]

---

Langkah 3: Hitung Nilai Ekspektasi \(\tilde{a}_j\)

Cari solusi persamaan kuadrat untuk setiap strata:

1. Strata <40 tahun:
   \[ (1-3)\tilde{a}_1^2 + [170 - 100 + 3(30 + 100)]\tilde{a}_1 - 100 \times 30 \times 3 = 0 \]
   \[ -2\tilde{a}_1^2 + 320\tilde{a}_1 - 9000 = 0 \]
   Solusi valid: \(\tilde{a}_1 = 30\) (dari rumus ABC).

2. Strata ≥40 tahun:
   \[ -2\tilde{a}_2^2 + 260\tilde{a}_2 - 21000 = 0 \]
   Solusi valid: \(\tilde{a}_2 = 70\).

---

Langkah 4: Hitung Varians

• Strata <40 tahun:
  \[ \tilde{b}_1 = 100-30=70, \tilde{c}_1 = 30-30=0, \tilde{d}_1 = 170-0=170 \]
  *Karena \(\tilde{c}_1=0\), varians tidak terdefinisi \(\Rightarrow\) Gunakan koreksi kontinuitas atau batasi nilai minimum 0.5.

(Misal: Set \(\tilde{c}_1=0.5\) untuk menghindari pembagian 0)

---

Langkah 5: Statistik Breslow-Day

\[ X_{\text{HBD}}^2 = \frac{(20-30)^2}{\text{Var}(a_1)} + \frac{(60-70)^2}{\text{Var}(a_2)} \]
(Nilai varians dihitung dengan koreksi, misal \(\text{Var}(a_1) \approx 15\))

Hasil:
- \(X_{\text{HBD}}^2 \approx 6.67\)
- \(df = 2-1 = 1\)
- Nilai kritis \(\chi^2_{1, 0.05} = 3.84\)

Keputusan:
- Karena \(6.67 > 3.84\), tolak \(H_0\) \(\Rightarrow\) OR tidak homogen antar strata.

---

Menghitung p-value

# Membuat array 3D (2 strata x 2x2)
data <- array(
  c(20, 60, 80, 40,  # Strata 1: <40 tahun (Ya/Tidak)
     10, 30, 90, 70), # Strata 2: ≥40 tahun (Ya/Tidak)
  dim = c(2, 2, 2),
  dimnames = list(
    Merokok = c("Ya", "Tidak"),
    Kanker = c("Ya", "Tidak"),
    Usia = c("<40 tahun", "≥40 tahun")
  )
)

# Install package jika belum ada
# install.packages("DescTools")

library(DescTools)

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

hasil_uji <- BreslowDayTest(data, correct = TRUE)
print(hasil_uji)

## 
##  Breslow-Day Test on Homogeneity of Odds Ratios (with Tarone correction)
## 
## data:  data
## X-squared = 0.7419, df = 1, p-value = 0.3891

Interpretasi Hasil Uji Breslow-Day

Output:
- Statistik Uji: \(X^2 = 0.7419\)
- Derajat Kebebasan: \(df = 1\)
- \(p\text{-value}\): \(0.3891\)

Kesimpulan:
- Karena \(p\text{-value} > 0.05\) (\(0.3891 \gg 0.05\)), gagal tolak \(H_0\).
- Artinya: Tidak ada bukti cukup untuk menyatakan bahwa odds ratio (OR) berbeda secara signifikan antar strata (\(\theta_{XY(1)} = \theta_{XY(2)}\)).

Implikasi:
1. Homogenitas Terjaga:
OR bersifat konsisten di semua strata variabel kontrol (\(Z\)).
Contoh: Jika \(Z\) adalah usia, efek merokok (\(X\)) pada kanker (\(Y\)) tidak berbeda antara kelompok usia.

2. Analisis Lanjutan:
   • Hasil ini mendukung penggunaan OR gabungan Mantel-Haenszel (\(\widehat{OR}_{MH}\)) untuk inferensi keseluruhan.
     
   • Tidak perlu melaporkan hasil per strata secara terpisah (kecuali ada alasan teoritis).

Catatan:
- Koreksi Tarone digunakan untuk mengurangi bias pada sampel kecil.
- Nilai \(X^2\) rendah (\(0.7419\)) menunjukkan perbedaan OR antar strata minimal.

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 utama:

1. Distribusi dari exponential family untuk variabel respons.

2. Fungsi link yang menghubungkan ekspektasi dari variabel respons ke kombinasi linear dari prediktor

3. Fungsi linear prediktor adalah \(\eta = X\boldsymbol{\beta}\).

Exponential Family

Distribusi termasuk dalam exponential family jika dapat ditulis dalam bentuk:

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

Contoh distribusi yang termasuk dalam exponential family: - Distribusi Normal
- Distribusi Binomial
- Distribusi Poisson
- Distribusi Gamma

Contoh Pembuktian: Distribusi Binomial

Fungsi probabilitas distribusi binomial adalah:

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

Kita tuliskan ulang dalam bentuk exponential family:

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

Dengan substitusi:

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

• \(b(\theta) = -n \log(1 - \pi)\)

• \(\phi = 1\)

Maka distribusi binomial termasuk dalam exponential family.

Model 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 dan 1, dengan menggunakan fungsi aktivasi sigmoid.

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

Regresi logistik memodelkan hubungan antara satu atau lebih variabel independen dan mengklasifikasikan data ke dalam kelas-kelas diskrit. Model ini banyak digunakan dalam pemodelan prediktif, di mana model memetakan probabilitas matematis apakah suatu catatan termasuk ke dalam kategori tertentu atau tidak.

Sebagai contoh, angka 0 dapat mewakili kelas negatif, dan angka 1 mewakili kelas positif. Regresi logistik sangat umum digunakan untuk masalah klasifikasi biner, di mana hasil hanya memiliki dua kemungkinan kategori (0 dan 1).

Beberapa contoh penerapan klasifikasi biner di mana respons biner diharapkan antara lain:

• Menentukan probabilitas serangan jantung : Dengan bantuan model logistik, tenaga medis dapat menentukan hubungan antara variabel-variabel seperti berat badan, olahraga, dan sebagainya, untuk memprediksi apakah seseorang berisiko mengalami serangan jantung atau komplikasi medis lainnya.

• Kemungkinan diterima di universitas: Sistem aplikasi dapat memperkirakan probabilitas seorang siswa untuk diterima di universitas tertentu atau program studi tertentu dengan menganalisis hubunganantara variabel-variabel penentu seperti skor GRE, GMAT, atau TOEFL

• Mengidentifkasi email spam: Kotak masuk email diflter untuk menentukan apakah suatu email merupakan komunikasi promosi atau spam dengan cara memahami variabel prediktor dan menerapkan algoritma regresi logistik untuk memeriksa keasliannya

---

Keunggulan Utama Regresi Logistik

1. Cocok untuk data yang dapat dipisahkan secara linear: Data yang dapat dipisahkan secara linear sangat sesuai untuk digunakan pada model ini.
2. Dapat digunakan untuk probabilitas: Regresi logistik tidak hanya memberikan klasifikasi, tetapi juga memberikan estimasi probabilitas.
3. Penafsiran koefisien yang jelas: Koefisien dari regresi logistik memberikan informasi mengenai hubungan antara prediktor dan peluang terjadinya kejadian.
4. Tidak mengasumsikan distribusi normal dari variabel independen, berbeda dengan regresi linear klasik.

---

Fungsi Aktivasi Sigmoid

Model regresi logistik menggunakan fungsi sigmoid sebagai fungsi aktivasi, yang memiliki bentuk:

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

Penjelasan klasifikasi berdasarkan fungsi sigmoid:

- Jika hasil fungsi sigmoid lebih dari 0.5, maka output dianggap sebagai 1 (kelas positif).

- Jika hasilnya kurang dari 0.5, maka output diklasifikasikan sebagai 0 (kelas negatif).

- Jika grafik menuju ke arah negatif ekstrem, maka nilai prediksi \(y\) akan menjadi 0, dan sebaliknya.

Dengan kata lain, jika keluaran fungsi sigmoid adalah 0.65, maka itu berarti terdapat peluang sebesar 65% bahwa peristiwa tersebut akan terjadi — misalnya dalam kasus pelemparan koin. Jika nilai prediksi lebih besar dari ambang batas (misalnya 0.5), maka diklasifikasikan sebagai 1. Jika kurang dari 0.5, diklasifikasikan sebagai 0.

Simulasi data regresi logistik :

set.seed(33)
n <- 50
x <- seq(-5, 5, length.out = n)
log_odds <- -1 + 2 * x
prob <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, prob)
data <- data.frame(x = x, y = y, prob = prob)

plot(x, y, pch = 16, col = "gray60",
     xlab = "X", ylab = "Y / Probabilitas",
     main = "Regresi Logistik (n = 50)")
lines(x, prob, col = "blue", lwd = 2)
abline(h = 0.5, col = "red", lty = 2)
legend("topleft",
       legend = c("Data Biner", "Kurva Logistik", "Ambang 0.5"),
       col = c("gray60", "blue", "red"),
       pch = c(16, NA, NA),
       lty = c(NA, 1, 2),
       lwd = c(NA, 2, 1),
       pt.cex = 1.5,
       bty = "n")

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

Spesifikasi Model

Fungsi link (logit): Fungsi link function logit dapat dinyatakan dalam bentuk berikut:

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

Model regresi logistik:

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

Fungsi inverse link:

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

---

Estimasi Parameter

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

Fungsi log-likelihood untuk regresi logistik:

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

Dengan:

\[ \pi_i = \frac{\exp(\boldsymbol{x}_i^\top \boldsymbol{\beta})}{1 + \exp(\boldsymbol{x}_i^\top \boldsymbol{\beta})} \]

Estimasi dilakukan melalui iterasi Newton-Raphson atau algoritma Fisher Scoring.

Simulasi menggunakan syntax R :

# Set seed untuk replikasi
set.seed(33)

# Jumlah data
n <- 75

# Variabel prediktor: x ~ distribusi normal
x <- rnorm(n)

# Probabilitas berdasarkan fungsi logit
p <- 1 / (1 + exp(-(-0.5 + 2 * x)))

# Respons biner: y ~ distribusi binomial
y <- rbinom(n, 1, p)

# Gabungkan jadi data frame
data <- data.frame(x, y)

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

# Tampilkan ringkasan model
summary(model)

## 
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3782     0.2705  -1.398 0.162145    
## x             1.4694     0.4098   3.586 0.000336 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 102.889  on 74  degrees of freedom
## Residual deviance:  83.007  on 73  degrees of freedom
## AIC: 87.007
## 
## Number of Fisher Scoring iterations: 4

plot(data$x, data$y, pch = 16, col = "grey")
curve(predict(model, newdata = data.frame(x = x), type = "response")[order(x)],
add = TRUE, col = "blue", lwd = 2)

Contoh 2 :

Regresi logistik digunakan untuk memodelkan hubungan antara satu atau lebih variabel independen dengan variabel dependen yang bersifat biner (0/1). Model ini menghasilkan probabilitas kejadian berdasarkan fungsi logit (sigmoid).

set.seed(33)
n <- 200
x <- rnorm(n)
log_odds <- -0.3 + 2 * x
prob <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, prob)
data <- data.frame(y = y, x = x)

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

## 
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.2826     0.1763  -1.603    0.109    
## x             1.7942     0.2746   6.533 6.45e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.54  on 199  degrees of freedom
## Residual deviance: 196.56  on 198  degrees of freedom
## AIC: 200.56
## 
## Number of Fisher Scoring iterations: 5

Interpretasi Koefsien :

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

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

• Nilai p: menunjukkan signifkansi statistik dari prediktor

Koefsien log-odds dapat ditransformasikan ke odds ratio :

exp(coef(model))

## (Intercept)           x 
##   0.7538288   6.0146777

Prediksi dan Visualisasi Kurva Lagi

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

Evaluasi Model :

Dapat melihat akurasi model berdasarkan klasifikasi threshold 0.5.

# Klasifikasi dan confusion matrix
data$pred_class <- ifelse(data$pred > 0.5, 1, 0)
table(Predicted = data$pred_class, Actual = data$y)

##          Actual
## Predicted  0  1
##         0 83 28
##         1 23 66

Kesimpulan

• Model regresi logistik mampu memprediksi probabilitas berdasarkan hubungan antara variabel x dan hasil biner y.

• Kurva logit menunjukkan bentuk S yang memisahkan dua kelas secara probabilistik.

• Koefsien dapat diinterpretasikan dalam bentuk log-odds atau odds ratio.

• Evaluasi model menunjukkan performa klasifkasi yang baik jika p-value signifkan dan prediksi sesuai dengan observasi.

Contoh 3

Dalam contoh ini, kita akan menggunakan data asli mtcars dari R base package, dengan mengubah variabel target menjadi biner. Kita akan memodelkan apakah mobil memiliki transmisi otomatis (0) atau manual (1) berdasarkan variabel prediktor seperti mpg (mil per galon) dan hp (tenaga kuda).

data(mtcars)
df <- mtcars
df$am <- factor(df$am, labels = c("Automatic", "Manual"))

summary(df[, c("am", "mpg", "hp")])

##          am          mpg              hp       
##  Automatic:19   Min.   :10.40   Min.   : 52.0  
##  Manual   :13   1st Qu.:15.43   1st Qu.: 96.5  
##                 Median :19.20   Median :123.0  
##                 Mean   :20.09   Mean   :146.7  
##                 3rd Qu.:22.80   3rd Qu.:180.0  
##                 Max.   :33.90   Max.   :335.0

model_real <- glm(am ~ mpg + hp, data = df, family = binomial)
summary(model_real)

## 
## Call:
## glm(formula = am ~ mpg + hp, family = binomial, data = df)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -33.60517   15.07672  -2.229   0.0258 *
## mpg           1.25961    0.56747   2.220   0.0264 *
## hp            0.05504    0.02692   2.045   0.0409 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 43.230  on 31  degrees of freedom
## Residual deviance: 19.233  on 29  degrees of freedom
## AIC: 25.233
## 
## Number of Fisher Scoring iterations: 7

Interpretasi Koefsien

• Intercept: log-odds mobil menggunakan transmisi manual saat mpg dan hp bernilai nol.

• mpg: setiap kenaikan 1 mil/galon meningkatkan peluang mobil bertransmisi manual.

• hp: setiap kenaikan 1 tenaga kuda menurunkan peluang mobil bertransmisi manual.

• p-value < 0.05 menunjukkan signifkansi statistik.

Odds ratio dapat dihitung sebagai :

exp(coef(model_real))

##  (Intercept)          mpg           hp 
## 2.543664e-15 3.524063e+00 1.056588e+00

Interpretasi:

• Odds ratio mpg > 1 → semakin irit, semakin besar peluang bertransmisi manual.

• Odds ratio hp < 1 → semakin besar tenaga, semakin kecil peluang bertransmisi manual. Visualisasi Prediksi.

df$prob <- predict(model_real, type = "response")
df$pred_class <- ifelse(df$prob > 0.5, "Manual", "Automatic")
library(ggplot2)
ggplot(df, aes(x = mpg, y = hp, color = pred_class, shape = am)) +
geom_point(size = 3) +
labs(title = "Prediksi Regresi Logistik: Transmisi Mobil",
x = "Mil per Galon (mpg)",
y = "Tenaga Kuda (hp)",
color = "Prediksi", shape = "Fakta") +
theme_minimal()

Evaluasi Model

# Confusion matrix
table(Predicted = df$pred_class, Actual = df$am)

##            Actual
## Predicted   Automatic Manual
##   Automatic        16      3
##   Manual            3     10

Interpretasi :

Baris dan kolom menunjukkan prediksi dan aktual.

Diagonal = prediksi benar, lainnya = kesalahan klasifkasi.

Kesimpulan :

• Model berhasil memprediksi jenis transmisi mobil berdasarkan mpg dan hp.

• Koefsien mpg positif dan signifkan → irit bahan bakar → kemungkinan lebih besar manual.

• Koefsien hp negatif → mobil lebih bertenaga cenderung otomatis.

• Model memiliki interpretasi yang kuat dan dapat divisualisasikan dengan baik.

Model 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.

Distribusi Poisson memiliki fungsi probabilitas:

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

Kita dapat menuliskan bentuk ini dalam format exponential family:

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

dengan:

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

Maka distribusi Poisson termasuk dalam exponential family.

Fungsi Link

Fungsi link kanonik untuk distribusi Poisson adalah fungsi logaritma:

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

Sehingga modelnya menjadi:

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

dan fungsi inverse link:

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

Estimasi Parameter

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

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

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

Contoh 1 :

# Simulasi data dengan hubungan yang kuat antara x dan y
set.seed(33)
n <- 200
x <- rnorm(n)
lambda <- exp(0.5 + 1.2 * x)  # hubungan kuat antara x dan log(mu)
y <- rpois(n, lambda)
data <- data.frame(y, x)

# Estimasi Regresi Poisson
poisson_model <- glm(y ~ x, data = data, family = poisson)
summary(poisson_model)

## 
## Call:
## glm(formula = y ~ x, family = poisson, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.54559    0.05640   9.673   <2e-16 ***
## x            1.19564    0.02834  42.193   <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: 1971.97  on 199  degrees of freedom
## Residual deviance:  230.57  on 198  degrees of freedom
## AIC: 680.13
## 
## Number of Fisher Scoring iterations: 5

plot(data$x, data$y, pch = 16, col = "darkgray", main = "Data dan Hasil Prediksi")
newdata <- data.frame(x = seq(min(x), max(x), length.out = 100))
pred <- predict(poisson_model, newdata = newdata, type = "response")
lines(newdata$x, pred, col = "blue", lwd = 2)

Diagnostic dan Overdispersion
Salah satu asumsi penting dari model Poisson adalah bahwa mean dan varians dari variabel response adalah sama \(E[Y] = Var[Y]\). Jika varians lebih besar dari mean, maka terjadi overdispersion.

Untuk mendeteksi overdispersion:

dispersion <- sum(residuals(poisson_model, type = "pearson")^2) / poisson_model$df.residual
dispersion

## [1] 0.9730813

Jika nilai dispersion > 1, maka overdispersion mungkin terjadi dan model alternatif seperti Negative Binomial Regression dapat digunakan.

Model regresi Poisson adalah alat penting untuk menganalisis data cacah. Ia memberikan hubungan log-linear antara prediktor dan rata-rata kejadian. Namun, perlu diperhatikan kemungkinan overdispersion dalam penerapannya.

Contoh 2:
Regresi Poisson digunakan untuk memodelkan data cacah (count data), yaitu data response berupa bilangan bulat non-negatif. Dalam contoh ini, kita akan menggunakan dataset warpbreaks dari R base package yang berisi jumlah patahan benang pada mesin tenun.

# Load data
data("InsectSprays")
head(InsectSprays)

summary(InsectSprays)

##      count       spray 
##  Min.   : 0.00   A:12  
##  1st Qu.: 3.00   B:12  
##  Median : 7.00   C:12  
##  Mean   : 9.50   D:12  
##  3rd Qu.:14.25   E:12  
##  Max.   :26.00   F:12

• Count : jumlah serangga tersisa setelah disemprot

• Spray : Jenis Semprotan

# Fit model Poisson
poisson_model <- glm(count ~ spray, family = poisson, data = InsectSprays)
summary(poisson_model)

## 
## Call:
## glm(formula = count ~ spray, family = poisson, data = InsectSprays)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.67415    0.07581  35.274  < 2e-16 ***
## sprayB       0.05588    0.10574   0.528    0.597    
## sprayC      -1.94018    0.21389  -9.071  < 2e-16 ***
## sprayD      -1.08152    0.15065  -7.179 7.03e-13 ***
## sprayE      -1.42139    0.17192  -8.268  < 2e-16 ***
## sprayF       0.13926    0.10367   1.343    0.179    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 409.041  on 71  degrees of freedom
## Residual deviance:  98.329  on 66  degrees of freedom
## AIC: 376.59
## 
## Number of Fisher Scoring iterations: 5

Interpretasi Koefisien

• \(\text{Intercept} = 2.67415\): nilai log dari rata-rata jumlah serangga untuk spray A (kategori referensi).

• \(\text{sprayB} = 0.05588\): rata-rata jumlah serangga untuk spray B sekitar \(e^{0.05588} \approx 1.06\) kali lebih banyak dibanding spray A (tidak signifikan, \(p = 0.597\)).

• \(\text{sprayC} = -1.90481\): rata-rata jumlah serangga untuk spray C sekitar \(e^{-1.90481} \approx 0.15\) kali dibanding spray A (sangat signifikan, \(p < 0.001\)).

• \(\text{sprayD} = -1.08152\): rata-rata jumlah serangga untuk spray D sekitar \(e^{-1.08152} \approx 0.34\) kali dibanding spray A (sangat signifikan, \(p < 0.001\)).

• \(\text{sprayE} = -1.42139\): rata-rata jumlah serangga untuk spray E sekitar \(e^{-1.42139} \approx 0.24\) kali dibanding spray A (sangat signifikan, \(p < 0.001\)).

• \(\text{sprayF} = 0.13926\): rata-rata jumlah serangga untuk spray F sekitar \(e^{0.13926} \approx 1.15\) kali dibanding spray A (tidak signifikan, \(p = 0.179\)).

# Lihat rate ratio (eksponen dari koefisien)
exp(coef(poisson_model))

## (Intercept)      sprayB      sprayC      sprayD      sprayE      sprayF 
##  14.5000000   1.0574713   0.1436782   0.3390805   0.2413793   1.1494253

• Jika OR < 1 → menurunkan jumlah serangga yang tersisa • Jika OR > 1 → meningkatkan jumlah serangga yang tersisa

# Tambahkan prediksi ke dalam data
InsectSprays$predicted <- predict(poisson_model, type = "response")
# Visualisasi
library(ggplot2)
ggplot(InsectSprays, aes(x = spray, y = count, color = spray)) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  geom_point(aes(y = predicted), shape = 18, size = 3, color = "black") +
  labs(
    title = "Prediksi Jumlah Serangga berdasarkan Jenis Semprotan",
    x = "Jenis Semprotan",
    y = "Jumlah Serangga Tersisa",
    color = "Spray"
  ) +
  theme_minimal()

Evaluasi Model

# Plot residual dari model Poisson untuk InsectSprays
plot(
  poisson_model$residuals,
  main = "Residual Plot",
  ylab = "Residual",
  xlab = "Index",
  pch = 19,
  col = "blue"
)

# Tambahkan garis horizontal di 0
abline(h = 0, col = "red", lty = 2)

Interpretasi :

• Regresi Poisson menunjukkan bahwa jenis semprotan berpengaruh terhadap jumlah serangga yang tersisa.

• Koefisien negatif menunjukkan penurunan rata-rata jumlah serangga.

• Model dapat digunakan untuk memprediksi jumlah serangga yang tersisa berdasarkan jenis semprotan yang digunakan.

• Beberapa jenis semprotan (C, D, dan E) menunjukkan efek yang signifikan dalam menurunkan jumlah serangga dibandingkan semprotan

Inferensi GLM

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

Ekspektasi dan Varians dalam GLM

1. Ekspektasi Estimator

Ekspektasi menunjukkan apakah suatu estimator tak bias, yaitu:

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

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

2. Varians Estimator

Varians menunjukkan presisi dari estimasi parameter:

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

di mana \(W\) adalah matriks bobot yang tergantung pada distribusi dan fungsi link.

\(\textit{Distribusi Asimptotik Estimator}\)

Dengan ukuran sampel besar:

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

Distribusi ini adalah dasar dari:

• \(\textit{Uji Wald}\)
• \(\textit{Confidence interval}\)
• \(\textit{P-value}\)

\(\textit{Varians dalam GLM Tidak Konstan}\)

Tidak seperti regresi linear (OLS) yang mengasumsikan homoskedastisitas:

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

Dalam GLM:

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

• \(\phi\) = parameter dispersi
• \(V(\mu)\) = fungsi varians

Contoh:

• \(\textit{Poisson}\): \(V(\mu) = \mu\)
• \(\textit{Binomial}\): \(V(\mu) = \mu(1 - \mu)\)

Contoh regressi poisson

# Simulasi data biner dengan prediktor tunggal
set.seed(123)
x <- rnorm(100)
logit_p <- -1 + 2.5*x  # Fungsi link logit
p <- plogis(logit_p)    # Transformasi ke probabilitas
y <- rbinom(100, 1, p)  # Data biner

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

## 
## Call:
## glm(formula = y ~ x, family = binomial)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.7479     0.3024  -2.473   0.0134 *  
## x             2.8406     0.5638   5.038  4.7e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 136.663  on 99  degrees of freedom
## Residual deviance:  77.592  on 98  degrees of freedom
## AIC: 81.592
## 
## Number of Fisher Scoring iterations: 6

Interpretasi:

• Koefisien \(\hat{\beta}_0 = -0.7479\) menunjukkan log-odds kejadian \(y=1\) ketika \(x=0\). Nilai \(p = 0.0134 < 0.05\) berarti intercept signifikan pada taraf 5%.

• Koefisien \(\hat{\beta}_1 = 2.8466\) berarti setiap kenaikan 1 unit pada \(x\) meningkatkan log-odds \(y=1\) sebesar 2.8466. Nilai \(p = 4.7 \times 10^{-7} < 0.05\) menunjukkan \(x\) sangat signifikan.

Penurunan Ekspektasi dan Varians dalam Generalized Linear Models

Konsep Dasar Ekspektasi

Berdasarkan pendekatan fungsi momen, nilai harapan dapat dinyatakan sebagai:

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

Untuk distribusi dari keluarga eksponensial, fungsi likelihood memiliki bentuk umum:

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

Atau dalam notasi alternatif:

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

Turunan pertama fungsi log-likelihood memberikan:

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

Dengan mengambil nilai harapannya diperoleh:

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

Sehingga hubungan fundamentalnya adalah:

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

Penurunan Varians

Melalui turunan kedua fungsi log-likelihood:

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

Diperoleh hubungan varians sebagai:

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

Metode Estimasi Parameter

Estimasi Kemungkinan Maksimum (MLE)
Proses estimasi parameter melalui optimasi fungsi likelihood:
\[ \ell(\theta) = \log L(\theta) \]
Langkah implementasi:
1. Cari titik stasioner dimana turunan pertama bernilai nol:
\[ U(\theta) = \frac{\partial\ell}{\partial\theta} = 0 \]
2. Verifikasi titik maksimum dengan uji turunan kedua:
\[ \frac{\partial^2\ell}{\partial\theta^2} < 0 \]

Algoritma Numerik
Untuk model non-linear, digunakan pendekatan iteratif:

Newton-Raphson:
\[ \beta^{(t+1)} = \beta^{(t)} - H^{-1}(\beta^{(t)})U(\beta^{(t)}) \]
dengan komponen:
- \(U(\beta)\): Vektor gradien
- \(H(\beta)\): Matriks Hessian

Fisher Scoring:
Modifikasi menggunakan matriks informasi Fisher \(I(\beta)\)

IRLS:
Solusi berbasis kuadrat terkecil terboboti:
\[ \beta^{(t+1)} = (X^TWX)^{-1}X^TWz \]

Diagnostic Model GLM

Ukuran Devians
Membandingkan model current vs saturated:
\[ D = 2\sum \left[y_i\log\left(\frac{y_i}{\hat{\mu}_i}\right) - (y_i - \hat{\mu}_i)\right] \]
Interpretasi:
- Devians kecil → model baik
- Devians besar → perlu revisi spesifikasi

Statistik Chi square Pearson
Mengukur goodness-of-fit:
\[ X^2 = \sum \frac{(y_i - \hat{\mu}_i)^2}{V(\hat{\mu}_i)} \]
Untuk data terkelompok mengikuti distribusi \(\chi^2\)

Analisis Residual
Diagnostik melalui:
- Residual deviance: \(r_i^D = \text{sign}(y_i - \hat{\mu}_i)\sqrt{d_i}\)
- Residual Pearson: \(r_i^P = \frac{y_i - \hat{\mu}_i}{\sqrt{V(\hat{\mu}_i)}}\)

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(X\beta)}{1+\exp(X\beta)} \]

Log-likelihood untuk n observasi :
\[ \ell(\beta) = \sum [y_i\log\pi_i + (1-y_i)\log(1-\pi_i)] \]

Estimasi dengan Newton-Raphson

Metode Newton-Raphson digunakan untuk mencari nilai parameter \(\beta\) yang memaksimalkan fungsi

loglikelihood pada model regresi logistik.

Prosedur Estimasi

1. Turunan Pertama (Score Function)
\[ U(\beta) = X^T(y - \pi) \]
2. Turunan Kedua (Hessian Matrix)
\[ W = \text{diag}(\pi_i(1-\pi_i)) \]
3. Iterasi Newton-Raphson
\[ \beta^{(t+1)} = \beta^{(t)} + (X^TWX)^{-1}X^T(y - \pi^{(t)}) \]

Estimasi MLE dengan Newton-Raphson

# Simulasi Data Poisson Regression
set.seed(33)
n <- 100  # Jumlah observasi
x <- rnorm(n, mean = 1, sd = 0.5)
beta_true <- c(0.5, 1.2)  # Intercept dan slope
X <- cbind(1, x)  # Matriks desain

# Membuat nilai prediksi
eta <- X %*% beta_true
lambda <- exp(eta)  # Fungsi link
y <- rpois(n, lambda)  # Data respon

# Implementasi Newton-Raphson Manual
beta <- c(0, 0)  # Nilai awal
tol <- 1e-6
max_iter <- 100
converged <- FALSE

for (i in 1:max_iter) {
  eta <- X %*% beta
  lambda_hat <- exp(eta)
  
  # Gradient (turunan pertama)
  U <- t(X) %*% (y - lambda_hat)
  
  # Hessian (turunan kedua)
  W <- diag(as.numeric(lambda_hat))
  H <- -t(X) %*% W %*% X
  
  # Update parameter
  beta_new <- beta - solve(H) %*% U
  
  # Cek konvergensi
  if (max(abs(beta_new - beta)) < tol) {
    converged <- TRUE
    break
  }
  
  beta <- beta_new
}
beta

##        [,1]
##   0.5467443
## x 1.2011478

Inferensi Parameter

1. Uji Wald

Tujuan Pengujian

Digunakan untuk menguji signifikansi statistik dari parameter individu \(\beta_j\) dalam model regresi logistik:

• Hipotesis nol (\(H_0\)): \(\beta_j = 0\)
  
• Hipotesis alternatif (\(H_1\)): \(\beta_j \neq 0\)

Landasan Teoretis

Berdasarkan sifat estimator MLE, \(\hat{\beta}_j\) mengikuti distribusi normal asimtotik:

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

Di bawah asumsi \(H_0\) benar (\(\beta_j = 0\)), statistik uji:

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

Statistik Wald diperoleh dengan mengkuadratkan \(Z\):

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

Langkah-langkah pengujian:

1. Hitung estimasi parameter \(\hat{\beta}_j\) melalui MLE
2. Estimasi standard error \(SE(\hat{\beta}_j)\) dari matriks informasi Fisher
3. Hitung statistik uji \(W\)
4. Bandingkan dengan nilai kritis \(\chi^2_1\) atau hitung p-value:
   • \(p = P(\chi^2_1 > W_{\text{obs}})\)
5. Tolak \(H_0\) jika \(p < \alpha\) (biasanya \(\alpha = 0.05\))

## Contoh Kasus Uji Wald: Pengaruh Usia dan Status Merokok terhadap Kanker Paru

set.seed(456)
n <- 200
usia <- round(rnorm(n, mean=50, sd=10))
merokok <- rbinom(n, 1, 0.4) # 1=perokok, 0=tidak
log_odds <- -3 + 0.05*usia + 1.8*merokok
prob <- plogis(log_odds) # inverse logit
kanker <- rbinom(n, 1, prob)
data_kanker <- data.frame(usia, merokok, kanker)

model_kanker <- glm(kanker ~ usia + merokok, 
                   data = data_kanker,
                   family = binomial)
summary(model_kanker)

## 
## Call:
## glm(formula = kanker ~ usia + merokok, family = binomial, data = data_kanker)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.72880    0.83925  -3.251  0.00115 ** 
## usia         0.04881    0.01612   3.028  0.00246 ** 
## merokok      1.45796    0.32235   4.523  6.1e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 272.74  on 199  degrees of freedom
## Residual deviance: 241.70  on 197  degrees of freedom
## AIC: 247.7
## 
## Number of Fisher Scoring iterations: 4

Interpretasi :

1. Intercept (\(\beta_0\) = -2.72880):
   • Log-odds kanker ketika usia=0 dan bukan perokok
   • Nilai negatif menunjukkan probabilitas dasar rendah
   • Signifikan secara statistik (p=0.00115)
2. Usia (\(\beta_1\) = 0.04881):
   • Setiap kenaikan 1 tahun usia meningkatkan log-odds kanker sebesar 0.04881
   • Odds Ratio: \(\exp(0.04881) \approx 1.050\) (kenaikan 5% risiko per tahun)
   • Sangat signifikan (p=0.00246)
3. Merokok (\(\beta_2\) = 1.45796):
   • Perokok memiliki log-odds kanker lebih tinggi 1.45796 dibanding non-perokok
   • Odds Ratio: \(\exp(1.45796) \approx 4.30\) (risiko 4.3 kali lebih tinggi)
   • Sangat signifikan (p=6.1×10⁻⁶)

beta_hat <- coef(model)["x"]
se_beta <- summary(model)$coefficients["x", "Std. Error"]
Z <- beta_hat / se_beta
Z

##       x 
## 6.53289

Wald_stat <- Z^2
Wald_stat

##        x 
## 42.67865

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

##            x 
## 6.451273e-11

Interpretasi

Jika p-value < 0.05, maka koefsien signifkan → variabel prediktor berpengaruh.

Jika p-value > 0.05, maka tidak ada cukup bukti untuk menolak H

Kesimpulan Uji Wald didasarkan pada rasio antara estimasi parameter dan standar error-nya. Dengan menaikkan nilai

Z menjadi kuadrat (Z²), kita memperoleh distribusi Chi-Square untuk pengujian hipotesis parameter

individual dalam model regresi logistik.

2. Uji Likelihood Ratio (Chi-Squdare)

set.seed(789)
n <- 500
pendidikan <- sample(1:4, n, replace=TRUE)  # 1=SD, 2=SMP, 3=SMA, 4=S1+
pendapatan <- rnorm(n, mean=50, sd=10)  # Pendapatan (juta/tahun)
log_odds <- -4 + 0.8*pendidikan + 0.05*pendapatan
prob <- plogis(log_odds)
kepemilikan_rumah <- rbinom(n, 1, prob)
data_rumah <- data.frame(pendidikan, pendapatan, kepemilikan_rumah)

# Model null (tanpa prediktor)
model_null <- glm(kepemilikan_rumah ~ 1, data=data_rumah, family=binomial)

# Model lengkap (dengan pendidikan + pendapatan)
model_full <- glm(kepemilikan_rumah ~ pendidikan + pendapatan, 
                 data=data_rumah, family=binomial)

# Likelihood Ratio Test
anova(model_null, model_full, test="Chisq")

Evaluasi kebaikan model :

AIC(model_full)

## [1] 574.1077

Semakin kecil nilai AIC semakin baik model

BIC(model_full)

## [1] 586.7516

AIC dan BIC digunakan untuk mengevaluasi kompleksitas dan kecocokan model.

Metode Estimasi dan Inferensi untuk Regresi Poisson

Konsep Dasar Model

Digunakan untuk data count (diskrit positif) dengan asumsi: - \(Y_i \sim \text{Poisson}(\lambda_i)\) - Fungsi massa probabilitas: \[ P(Y_i=y_i) = \frac{e^{-\lambda_i}\lambda_i^{y_i}}{y_i!} \]

Spesifikasi Model

Hubungan antara prediktor dan respons melalui fungsi link log: \[ \log(\lambda_i) = \mathbf{x}_i^\top\beta \quad \Rightarrow \quad \lambda_i = \exp(\mathbf{x}_i^\top\beta) \]

Estimasi Parameter

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

Algoritma IRLS (Iteratively Reweighted Least Squares): 1. Hitung linear predictor \(\eta_i = \mathbf{x}_i^\top\beta\) 2. Hitung nilai fitted \(\lambda_i = \exp(\eta_i)\) 3. Update parameter: \[ \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)\) (matriks bobot)

• \(\mathbf{z} = \eta + \frac{y - \lambda}{\lambda}\) (variabel respon adjusted)

Implementasi Manual dalam R

### Contoh Kasus: Model Frekuensi Kecelakaan Lalu Lintas

#### Simulasi Data
#Faktor prediktor:
# `jam_sibuk`: 1=jam sibuk (18.00-20.00), 0=lainnya
# `curah_hujan`: Tingkat curah hujan (mm)
# `kondisi_jalan`: 1=rusak, 0=baik

# Simulasi Data Kecelakaan
set.seed(123)
n <- 300
jam_sibuk <- rbinom(n, 1, 0.4)
curah_hujan <- round(runif(n, min=0, max=15), 1)
kondisi_jalan <- rbinom(n, 1, 0.3)

# True parameters
beta_true <- c(1.5, 0.8, 0.02, 0.5)

# Linear predictor
X <- cbind(1, jam_sibuk, curah_hujan, kondisi_jalan)
lambda <- exp(X %*% beta_true)

# Generate response
kecelakaan <- rpois(n, lambda)
data_accident <- data.frame(jam_sibuk, curah_hujan, kondisi_jalan, kecelakaan)

summary(data_accident)

##    jam_sibuk     curah_hujan     kondisi_jalan   kecelakaan    
##  Min.   :0.00   Min.   : 0.000   Min.   :0.0   Min.   : 0.000  
##  1st Qu.:0.00   1st Qu.: 4.025   1st Qu.:0.0   1st Qu.: 5.000  
##  Median :0.00   Median : 7.300   Median :0.0   Median : 8.000  
##  Mean   :0.39   Mean   : 7.481   Mean   :0.3   Mean   : 9.107  
##  3rd Qu.:1.00   3rd Qu.:11.400   3rd Qu.:1.0   3rd Qu.:12.000  
##  Max.   :1.00   Max.   :15.000   Max.   :1.0   Max.   :29.000

hist(data_accident$kecelakaan, main="Distribusi Frekuensi Kecelakaan", xlab="Jumlah Kecelakaan")

model_pois <- glm(kecelakaan ~ jam_sibuk + curah_hujan + kondisi_jalan,
                 family=poisson, data=data_accident)
summary(model_pois)

## 
## Call:
## glm(formula = kecelakaan ~ jam_sibuk + curah_hujan + kondisi_jalan, 
##     family = poisson, data = data_accident)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.438898   0.049696  28.954  < 2e-16 ***
## jam_sibuk     0.831283   0.039152  21.232  < 2e-16 ***
## curah_hujan   0.022440   0.004371   5.133 2.84e-07 ***
## kondisi_jalan 0.534102   0.038740  13.787  < 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: 938.41  on 299  degrees of freedom
## Residual deviance: 285.98  on 296  degrees of freedom
## AIC: 1458
## 
## Number of Fisher Scoring iterations: 4

# 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: 6.53289 
## Chi-Square: 42.67865 
## p-value: 6.451273e-11

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

## Warning in anova.glmlist(c(list(object), dotargs), dispersion = dispersion, :
## models with response '"y"' removed because response differs from model 1

Evaluasi Model dengan AIC dan BIC

AIC(model_null1)

## [1] 2419.529

BIC(model_null1)

## [1] 2422.828

Interpretasi Singkat Hasil Analisis:

1. Model Regresi Poisson:

   • Jam sibuk meningkatkan risiko kecelakaan 2.3× (p < 0.001).

   • Curah hujan (per 1mm) tingkatkan risiko 2.2% (p = 2.8e-07).

   • Jalan rusak tingkatkan risiko 71% (p < 0.001).

   • Kesesuaian model: Deviance turun drastis (938 → 286), menunjukkan prediktor sangat berpengaruh.

2. Likelihood Ratio Test:

   • Model dengan prediktor secara signifikan lebih baik daripada model null (p = 1.8e-11).

   • Penurunan deviance 1065.2 mengonfirmasi kekuatan prediktor.

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 kemahiran teknologi (Beginner < Intermediate < Advanced < Expert). Dalam analisis, variabel ordinal bisa diperlakukan :

  – Sebagai nominal dengan membuat dummy untuk setiap kategori. – Atau diperlakukan seperti rasio dengan memberi nilai tertentu (misalnya beginner = 1, intermediate = 2, advanced = 3, expert = 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

# Load the necessary library
library(tibble)

# 1. Set simulation parameters
n <- 333
set.seed(33)

# 2. Simulate independent variables
gender <- sample(c("Male", "Female"), n, replace = TRUE)
tech_savviness <- sample(
  c("Beginner", "Intermediate", "Advanced", "Expert"),
  n,
  replace = TRUE,
  prob = c(0.20, 0.45, 0.25, 0.10) # Assigning probabilities to each level
)

income <- rnorm(n, mean = 50, sd = 15)

# 3. Model the probability of success ('p')
logit_p <- -2.5 +
  0.5 * (gender == "Female") +
  0.7 * as.numeric(factor(tech_savviness, ordered = TRUE)) + # Higher level -> higher success
  0.03 * income

p <- 1 / (1 + exp(-logit_p))

# 4. Simulate the final 'success' outcome
success <- rbinom(n, 1, p)

# 5. Combine all variables into a single data frame (tibble)
sim_data <- tibble(success, gender, tech_savviness, income)

# 6. Display the first few rows of the new simulated data
head(sim_data)

Interpretasi : Dataset berisi 333 observasi dengan variabel gender, kemahiran teknologi, income, dan status keberhasilan.

Eksplorasi Data

sim_data %>%
dplyr::group_by(success) %>%
dplyr::summarise(
  Jumlah = dplyr::n(),
  Rata2_Income = mean(income)
)

Interpretasi : Distribusi jumlah sukses dan tidak sukses, serta rata-rata pendapatan pada masing-masing grup.

Perlakuan Variabel Ordinal

Treat Sebagai Nominal (Dummy)

library(dplyr)

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

## 
## Attaching package: 'dplyr'

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

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

sim_data_nominal <- sim_data %>%
mutate(
tech_savviness = factor(tech_savviness, levels = c("Beginner", "Intermediate", "Advanced", "Expert"))
)
model_nominal <- glm(success ~ gender + tech_savviness + income, data = sim_data, family = binomial)

# Menampilkan ringkasan model
summary(model_nominal)

## 
## Call:
## glm(formula = success ~ gender + tech_savviness + income, family = binomial, 
##     data = sim_data)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -1.93421    0.57942  -3.338 0.000843 ***
## genderMale                 -0.06008    0.26950  -0.223 0.823591    
## tech_savvinessBeginner      0.84971    0.35339   2.404 0.016195 *  
## tech_savvinessExpert        1.45863    0.47432   3.075 0.002104 ** 
## tech_savvinessIntermediate  2.10844    0.34308   6.146 7.96e-10 ***
## income                      0.03770    0.01016   3.709 0.000208 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 392.54  on 332  degrees of freedom
## Residual deviance: 337.14  on 327  degrees of freedom
## AIC: 349.14
## 
## Number of Fisher Scoring iterations: 4

Interpretasi Koefisien

Intercept (-1.93421)

• Ini adalah nilai log-odds dasar untuk individu yang menjadi kategori referensi di semua variabel: gender (kemungkinan Female), tech_savviness (level referensi yang tidak disebutkan, misal “Novice”), dan memiliki income = 0.

• Sangat signifikan (p=0.000843), artinya baseline (kelompok referensi) ini memiliki peluang sukses yang secara signifikan berbeda dari 50%.

• Odds Ratio = 0.14: Peluang sukses untuk kelompok referensi ini hanya 0.14 kali dari peluang gagalnya.

genderMale (-0.06008)

• Individu Male memiliki log-odds sukses 0.06008 lebih rendah dibandingkan individu Female (kategori referensi).

• Tidak signifikan (p=0.823591), artinya tidak terdapat perbedaan statistik yang nyata dalam peluang sukses antara laki-laki dan perempuan berdasarkan model ini.

• Odds Ratio = 0.94: Peluang sukses laki-laki adalah 94% dari peluang sukses perempuan. Namun, perbedaan ini tidak signifikan secara statistik.

tech_savvinessBeginner (0.84971)

• Individu dengan tingkat tech_savviness Beginner memiliki log-odds sukses 0.84971 lebih tinggi dibandingkan level referensi.

• Signifikan pada level 5% (p=0.016195), menunjukkan bahwa tingkat Beginner memiliki peluang sukses yang lebih tinggi secara signifikan dibandingkan level referensi.

• Odds Ratio = 2.34: Peluang sukses individu dengan tech savviness Beginner adalah 2.34 kali lebih besar daripada individu di level referensi.

tech_savvinessExpert (1.45863)

• Individu dengan tingkat tech_savviness Expert memiliki log-odds sukses 1.45863 lebih tinggi dibandingkan level referensi.

• Sangat signifikan (p=0.002104), menunjukkan bahwa tingkat Expert memiliki peluang sukses yang jauh lebih tinggi secara signifikan dibandingkan level referensi.

• Odds Ratio = 4.30: Peluang sukses individu dengan tech savviness Expert adalah 4.3 kali lebih besar daripada individu di level referensi.

tech_savvinessIntermediate (2.10844)

• Individu dengan tingkat tech_savviness Intermediate memiliki log-odds sukses 2.10844 lebih tinggi dibandingkan level referensi.

• Sangat signifikan (p=7.96e−10), menunjukkan bahwa tingkat Intermediate memiliki peluang sukses yang paling tinggi secara signifikan dibandingkan semua level lainnya.

• Odds Ratio = 8.24: Peluang sukses individu dengan tech savviness Intermediate adalah 8.24 kali lebih besar daripada individu di level referensi.

income (0.03770)

• Setiap kenaikan satu unit pada income akan meningkatkan log-odds sukses sebesar 0.03770.

• Sangat signifikan (p=0.000208), menunjukkan bahwa pendapatan memiliki pengaruh positif dan signifikan terhadap peluang sukses.

• Odds Ratio = 1.04: Untuk setiap kenaikan satu unit income, peluang untuk sukses meningkat sebesar 4%.

Interpretasi Goodness-of-Fit

• Null deviance (392.54): Ini adalah deviance dari model yang hanya menyertakan intercept (model tanpa prediktor sama sekali). Angka ini berfungsi sebagai baseline untuk mengukur performa model.

• Residual deviance (337.14): Ini adalah deviance dari model setelah semua variabel prediktor dimasukkan.

  • Adanya penurunan dari null deviance (392.54) ke residual deviance (337.14) menunjukkan bahwa variabel-variabel prediktor yang ada di dalam model secara kolektif berhasil memberikan informasi yang baik dan signifikan untuk menjelaskan variabel dependen.

• AIC (349.14):

  • AIC (Akaike Information Criterion) digunakan untuk membandingkan kualitas relatif dari beberapa model. Semakin kecil nilai AIC, semakin baik model tersebut dalam menyeimbangkan antara kecocokan data (goodness-of-fit) dan kompleksitas model (jumlah prediktor). Nilai ini akan sangat berguna jika Anda ingin membandingkan model ini dengan versi model yang lain.

Signifikansi Model

• Variabel income dan tech savviness (Intermediate, Advanced, Expert) signifikan memengaruhi peluang sukses.

• Variabel gender menunjukkan tren, dan signifikan pada taraf 5%.

• Tech savviness pada tingkat Beginner justru diasosiasikan dengan penurunan peluang sukses dibandingkan Beginner dalam data ini, yang mungkin terjadi akibat struktur data simulasi.

Kesimpulan Praktis

• Income dan tingkat tech savviness yang lebih tinggi (Intermediate, Advanced, Expert) merupakan prediktor kuat untuk peluang sukses.

• Gender berpotensi berpengaruh tetapi perlu data tambahan untuk konfirmasi.

• Model cukup baik dalam memprediksi dibandingkan model null, namun ada ruang untuk peninjauan dengan mungkin menambahkan variabel lain.

Treat Sebagai Rasio (Numeric Rank)

sim_data_numeric <- sim_data %>%
  mutate(
    tech_savviness_numeric = case_when(
      tech_savviness == "Beginner" ~ 1,
      tech_savviness == "Intermediate" ~ 2,
      tech_savviness == "Advanced" ~ 3,
      tech_savviness == "Expert" ~ 4
    )
  )

model_numeric <- glm(success ~ gender + tech_savviness_numeric + income, data = sim_data_numeric, family = binomial)

summary(model_numeric)

## 
## Call:
## glm(formula = success ~ gender + tech_savviness_numeric + income, 
##     family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)   
## (Intercept)             0.146760   0.568490   0.258   0.7963   
## genderMale             -0.075699   0.251707  -0.301   0.7636   
## tech_savviness_numeric -0.256533   0.138778  -1.849   0.0645 . 
## income                  0.029506   0.009237   3.194   0.0014 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 392.54  on 332  degrees of freedom
## Residual deviance: 378.51  on 329  degrees of freedom
## AIC: 386.51
## 
## Number of Fisher Scoring iterations: 4

Interpretasi Koefisien

Intercept (0.146760)

• Ini adalah log-odds dasar untuk individu dengan nilai referensi (misalnya, gender non-Male, tech_savviness_numeric = 0, dan income = 0).

• Tidak signifikan (p = 0.7963), artinya baseline ini tidak berbeda secara signifikan dari probabilitas 50%.

• Odds Ratio = 1.158: Peluang sukses dasar sedikit lebih tinggi dari baseline, namun tidak signifikan.

genderMale (-0.075699)

• Individu Male memiliki log-odds sukses yang lebih rendah sebesar 0.076 dibandingkan dengan referensi (misalnya, Female atau non-Male).

• p = 0.7636 (tidak signifikan), menunjukkan bahwa perbedaan peluang sukses antara laki-laki dan referensi tidak signifikan secara statistik.

• Odds Ratio = 0.927: Peluang sukses laki-laki sekitar 93% dari peluang sukses referensi, tetapi perbedaan ini tidak signifikan.

tech_savviness_numeric (-0.256533)

• Setiap kenaikan satu unit dalam tech_savviness_numeric mengurangi log-odds sukses sebesar 0.257.

• p = 0.0645 (mendekati signifikan pada taraf 10%), menunjukkan tren bahwa semakin tinggi tech_savviness, peluang sukses mungkin sedikit menurun, tetapi tidak signifikan pada taraf 5%.

• Odds Ratio = 0.774: Setiap kenaikan satu unit tech_savviness mengurangi peluang sukses sekitar 23%, tetapi hasil ini tidak signifikan secara kuat.

income (+0.029506)

• Setiap kenaikan 1 unit income (misal, 1 juta), meningkatkan log-odds sukses sebesar 0.030.

• p = 0.0014 (sangat signifikan), artinya pendapatan berhubungan positif dengan peluang sukses.

• Odds Ratio = 1.030: Setiap tambahan 1 juta meningkatkan peluang sukses sekitar 3.0%.

Interpretasi Goodness-of-Fit

• Null deviance (392.54): Deviance model tanpa prediktor (hanya intercept).

• Residual deviance (378.51): Deviance model setelah memasukkan prediktor (genderMale, tech_savviness_numeric, dan income).

  • Penurunan dari null deviance ke residual deviance menunjukkan bahwa prediktor memberikan tambahan informasi yang berarti dalam menjelaskan variasi data.

• AIC (386.51):

  • Nilai AIC yang lebih rendah dibandingkan model tanpa prediktor menunjukkan bahwa model ini memiliki keseimbangan yang lebih baik antara goodness-of-fit dan kompleksitas.

Signifikansi Model

• income (p = 0.0014) sangat signifikan dan memiliki pengaruh positif terhadap peluang sukses.

• tech_savviness_numeric (p = 0.0645) mendekati signifikan (p < 0.1) tetapi tidak mencapai signifikansi pada taraf 5%, menunjukkan tren negatif yang lemah.

• genderMale (p = 0.7636) tidak signifikan, artinya tidak ada bukti statistik yang kuat bahwa gender memengaruhi peluang sukses dalam model ini.

Kesimpulan Praktis

1. Pendapatan (income) memiliki dampak positif yang kuat: setiap kenaikan pendapatan meningkatkan peluang sukses secara signifikan.

2. Keterampilan teknologi (tech_savviness_numeric) menunjukkan tren bahwa semakin tinggi skor, peluang sukses sedikit menurun, tetapi efek ini tidak kuat secara statistik.

3. Gender (Male vs. referensi) tidak memberikan pengaruh yang signifikan terhadap peluang sukses dalam model ini.

4. Kinerja Model:

   • Model menunjukkan perbaikan dibanding model tanpa prediktor (null model).

   • Namun, masih ada ruang untuk peningkatan dengan menambahkan variabel lain yang mungkin lebih relevan atau memperbaiki pengukuran variabel yang ada.

Perbandingan model

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

## $AIC_Nominal
## [1] 349.1385
## 
## $AIC_Numeric
## [1] 386.5096

Kesimpulan & Rekomendasi

• Gunakan model nominal (AIC = 349.14) karena lebih baik dalam menjelaskan data.

• Hindari mengubah variabel kategorikal menjadi numerik jika tidak ada justifikasi teoretis, karena dapat mengurangi akurasi model.

Goodness of fit model

nullmod <- glm(success ~ 1, data = sim_data, family = binomial)
r2_nominal <- 1 - (logLik(model_nominal)/logLik(nullmod))
r2_numeric <- 1 - (logLik(model_numeric)/logLik(nullmod))
list(
McFadden_R2_Nominal = r2_nominal,
McFadden_R2_Numeric = r2_numeric
)

## $McFadden_R2_Nominal
## 'log Lik.' 0.1411395 (df=6)
## 
## $McFadden_R2_Numeric
## 'log Lik.' 0.03574678 (df=4)

library(ggplot2)  
sim_data_nominal <- sim_data_nominal %>% mutate(predicted = predict(model_nominal, type = "response"))
sim_data_numeric <- sim_data_numeric %>% mutate(predicted = predict(model_numeric, type = "response"))
# Plot untuk model nominal
sim_data_nominal %>%
ggplot(aes(x = income, y = predicted, color = tech_savviness)) +
geom_point(alpha = 0.6) +
labs(title = "Prediksi Probabilitas (Ordinal sebagai Nominal)", x = "Income", y = "Prediksi Probabilitas")

theme_minimal()

## List of 136
##  $ line                            :List of 6
##   ..$ colour       : chr "black"
##   ..$ linewidth    : num 0.5
##   ..$ linetype     : num 1
##   ..$ lineend      : chr "butt"
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ rect                            :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : chr "black"
##   ..$ linewidth    : num 0.5
##   ..$ linetype     : num 1
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ text                            :List of 11
##   ..$ family       : chr ""
##   ..$ face         : chr "plain"
##   ..$ colour       : chr "black"
##   ..$ size         : num 11
##   ..$ hjust        : num 0.5
##   ..$ vjust        : num 0.5
##   ..$ angle        : num 0
##   ..$ lineheight   : num 0.9
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ title                           : NULL
##  $ aspect.ratio                    : NULL
##  $ axis.title                      : NULL
##  $ axis.title.x                    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.75points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x.top                :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 2.75points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x.bottom             : NULL
##  $ axis.title.y                    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.75points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.y.left               : NULL
##  $ axis.title.y.right              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : num -90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 2.75points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text                       :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : chr "grey30"
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x                     :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.2points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x.top                 :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 2.2points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x.bottom              : NULL
##  $ axis.text.y                     :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 1
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.2points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.y.left                : NULL
##  $ axis.text.y.right               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 2.2points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.theta                 : NULL
##  $ axis.text.r                     :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0.5
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.2points 0points 2.2points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.ticks                      : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.ticks.x                    : NULL
##  $ axis.ticks.x.top                : NULL
##  $ axis.ticks.x.bottom             : NULL
##  $ axis.ticks.y                    : NULL
##  $ axis.ticks.y.left               : NULL
##  $ axis.ticks.y.right              : NULL
##  $ axis.ticks.theta                : NULL
##  $ axis.ticks.r                    : NULL
##  $ axis.minor.ticks.x.top          : NULL
##  $ axis.minor.ticks.x.bottom       : NULL
##  $ axis.minor.ticks.y.left         : NULL
##  $ axis.minor.ticks.y.right        : NULL
##  $ axis.minor.ticks.theta          : NULL
##  $ axis.minor.ticks.r              : NULL
##  $ axis.ticks.length               : 'simpleUnit' num 2.75points
##   ..- attr(*, "unit")= int 8
##  $ axis.ticks.length.x             : NULL
##  $ axis.ticks.length.x.top         : NULL
##  $ axis.ticks.length.x.bottom      : NULL
##  $ axis.ticks.length.y             : NULL
##  $ axis.ticks.length.y.left        : NULL
##  $ axis.ticks.length.y.right       : NULL
##  $ axis.ticks.length.theta         : NULL
##  $ axis.ticks.length.r             : NULL
##  $ axis.minor.ticks.length         : 'rel' num 0.75
##  $ axis.minor.ticks.length.x       : NULL
##  $ axis.minor.ticks.length.x.top   : NULL
##  $ axis.minor.ticks.length.x.bottom: NULL
##  $ axis.minor.ticks.length.y       : NULL
##  $ axis.minor.ticks.length.y.left  : NULL
##  $ axis.minor.ticks.length.y.right : NULL
##  $ axis.minor.ticks.length.theta   : NULL
##  $ axis.minor.ticks.length.r       : NULL
##  $ axis.line                       : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.line.x                     : NULL
##  $ axis.line.x.top                 : NULL
##  $ axis.line.x.bottom              : NULL
##  $ axis.line.y                     : NULL
##  $ axis.line.y.left                : NULL
##  $ axis.line.y.right               : NULL
##  $ axis.line.theta                 : NULL
##  $ axis.line.r                     : NULL
##  $ legend.background               : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ legend.margin                   : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
##   ..- attr(*, "unit")= int 8
##  $ legend.spacing                  : 'simpleUnit' num 11points
##   ..- attr(*, "unit")= int 8
##  $ legend.spacing.x                : NULL
##  $ legend.spacing.y                : NULL
##  $ legend.key                      : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ legend.key.size                 : 'simpleUnit' num 1.2lines
##   ..- attr(*, "unit")= int 3
##  $ legend.key.height               : NULL
##  $ legend.key.width                : NULL
##  $ legend.key.spacing              : 'simpleUnit' num 5.5points
##   ..- attr(*, "unit")= int 8
##  $ legend.key.spacing.x            : NULL
##  $ legend.key.spacing.y            : NULL
##  $ legend.frame                    : NULL
##  $ legend.ticks                    : NULL
##  $ legend.ticks.length             : 'rel' num 0.2
##  $ legend.axis.line                : NULL
##  $ legend.text                     :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.text.position            : NULL
##  $ legend.title                    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.title.position           : NULL
##  $ legend.position                 : chr "right"
##  $ legend.position.inside          : NULL
##  $ legend.direction                : NULL
##  $ legend.byrow                    : NULL
##  $ legend.justification            : chr "center"
##  $ legend.justification.top        : NULL
##  $ legend.justification.bottom     : NULL
##  $ legend.justification.left       : NULL
##  $ legend.justification.right      : NULL
##  $ legend.justification.inside     : NULL
##  $ legend.location                 : NULL
##  $ legend.box                      : NULL
##  $ legend.box.just                 : NULL
##  $ legend.box.margin               : 'margin' num [1:4] 0cm 0cm 0cm 0cm
##   ..- attr(*, "unit")= int 1
##  $ legend.box.background           : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ legend.box.spacing              : 'simpleUnit' num 11points
##   ..- attr(*, "unit")= int 8
##   [list output truncated]
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi TRUE
##  - attr(*, "validate")= logi TRUE

# Plot untuk model numeric
sim_data_numeric %>%
ggplot(aes(x = income, y = predicted, color = as.factor(tech_savviness_numeric))) +
geom_point(alpha = 0.6) +
labs(title = "Prediksi Probabilitas (Ordinal sebagai Numeric)", x = "Income", y = "Prediksi Probabilitas")

theme_minimal()

## List of 136
##  $ line                            :List of 6
##   ..$ colour       : chr "black"
##   ..$ linewidth    : num 0.5
##   ..$ linetype     : num 1
##   ..$ lineend      : chr "butt"
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ rect                            :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : chr "black"
##   ..$ linewidth    : num 0.5
##   ..$ linetype     : num 1
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ text                            :List of 11
##   ..$ family       : chr ""
##   ..$ face         : chr "plain"
##   ..$ colour       : chr "black"
##   ..$ size         : num 11
##   ..$ hjust        : num 0.5
##   ..$ vjust        : num 0.5
##   ..$ angle        : num 0
##   ..$ lineheight   : num 0.9
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ title                           : NULL
##  $ aspect.ratio                    : NULL
##  $ axis.title                      : NULL
##  $ axis.title.x                    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.75points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x.top                :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 2.75points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x.bottom             : NULL
##  $ axis.title.y                    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.75points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.y.left               : NULL
##  $ axis.title.y.right              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : num -90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 2.75points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text                       :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : chr "grey30"
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x                     :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.2points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x.top                 :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 2.2points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x.bottom              : NULL
##  $ axis.text.y                     :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 1
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.2points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.y.left                : NULL
##  $ axis.text.y.right               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 2.2points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.theta                 : NULL
##  $ axis.text.r                     :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0.5
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.2points 0points 2.2points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.ticks                      : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.ticks.x                    : NULL
##  $ axis.ticks.x.top                : NULL
##  $ axis.ticks.x.bottom             : NULL
##  $ axis.ticks.y                    : NULL
##  $ axis.ticks.y.left               : NULL
##  $ axis.ticks.y.right              : NULL
##  $ axis.ticks.theta                : NULL
##  $ axis.ticks.r                    : NULL
##  $ axis.minor.ticks.x.top          : NULL
##  $ axis.minor.ticks.x.bottom       : NULL
##  $ axis.minor.ticks.y.left         : NULL
##  $ axis.minor.ticks.y.right        : NULL
##  $ axis.minor.ticks.theta          : NULL
##  $ axis.minor.ticks.r              : NULL
##  $ axis.ticks.length               : 'simpleUnit' num 2.75points
##   ..- attr(*, "unit")= int 8
##  $ axis.ticks.length.x             : NULL
##  $ axis.ticks.length.x.top         : NULL
##  $ axis.ticks.length.x.bottom      : NULL
##  $ axis.ticks.length.y             : NULL
##  $ axis.ticks.length.y.left        : NULL
##  $ axis.ticks.length.y.right       : NULL
##  $ axis.ticks.length.theta         : NULL
##  $ axis.ticks.length.r             : NULL
##  $ axis.minor.ticks.length         : 'rel' num 0.75
##  $ axis.minor.ticks.length.x       : NULL
##  $ axis.minor.ticks.length.x.top   : NULL
##  $ axis.minor.ticks.length.x.bottom: NULL
##  $ axis.minor.ticks.length.y       : NULL
##  $ axis.minor.ticks.length.y.left  : NULL
##  $ axis.minor.ticks.length.y.right : NULL
##  $ axis.minor.ticks.length.theta   : NULL
##  $ axis.minor.ticks.length.r       : NULL
##  $ axis.line                       : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.line.x                     : NULL
##  $ axis.line.x.top                 : NULL
##  $ axis.line.x.bottom              : NULL
##  $ axis.line.y                     : NULL
##  $ axis.line.y.left                : NULL
##  $ axis.line.y.right               : NULL
##  $ axis.line.theta                 : NULL
##  $ axis.line.r                     : NULL
##  $ legend.background               : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ legend.margin                   : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
##   ..- attr(*, "unit")= int 8
##  $ legend.spacing                  : 'simpleUnit' num 11points
##   ..- attr(*, "unit")= int 8
##  $ legend.spacing.x                : NULL
##  $ legend.spacing.y                : NULL
##  $ legend.key                      : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ legend.key.size                 : 'simpleUnit' num 1.2lines
##   ..- attr(*, "unit")= int 3
##  $ legend.key.height               : NULL
##  $ legend.key.width                : NULL
##  $ legend.key.spacing              : 'simpleUnit' num 5.5points
##   ..- attr(*, "unit")= int 8
##  $ legend.key.spacing.x            : NULL
##  $ legend.key.spacing.y            : NULL
##  $ legend.frame                    : NULL
##  $ legend.ticks                    : NULL
##  $ legend.ticks.length             : 'rel' num 0.2
##  $ legend.axis.line                : NULL
##  $ legend.text                     :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.text.position            : NULL
##  $ legend.title                    :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.title.position           : NULL
##  $ legend.position                 : chr "right"
##  $ legend.position.inside          : NULL
##  $ legend.direction                : NULL
##  $ legend.byrow                    : NULL
##  $ legend.justification            : chr "center"
##  $ legend.justification.top        : NULL
##  $ legend.justification.bottom     : NULL
##  $ legend.justification.left       : NULL
##  $ legend.justification.right      : NULL
##  $ legend.justification.inside     : NULL
##  $ legend.location                 : NULL
##  $ legend.box                      : NULL
##  $ legend.box.just                 : NULL
##  $ legend.box.margin               : 'margin' num [1:4] 0cm 0cm 0cm 0cm
##   ..- attr(*, "unit")= int 1
##  $ legend.box.background           : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ legend.box.spacing              : 'simpleUnit' num 11points
##   ..- attr(*, "unit")= int 8
##   [list output truncated]
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi TRUE
##  - attr(*, "validate")= logi TRUE

Ringkasan koefisien model

#=======Ringkasan model nominal=========#
# Load required packages
library(broom)

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

library(knitr)

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

library(kableExtra)

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

## 
## Attaching package: 'kableExtra'

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

# Assuming your model is named 'model_nominal'
# Example model (replace with your actual model):
# model_nominal <- glm(formula, data = your_data, family = binomial)

# Create the formatted table
tidy(model_nominal) %>%
  kable(
    format = "html",
    caption = "Ringkasan Koefisien Model dengan Ordinal sebagai Nominal",
    align = c("l", "r", "r", "r", "r"),
    col.names = c("Term", "Estimate", "Std. Error", "Statistic", "p-value"),
    digits = c(NA, 6, 6, 4, 7)
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE,
    position = "left"
  ) %>%
  add_header_above(c(" " = 1, "Coefficients" = 4)) %>%
  footnote(
    general = "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
    general_title = "Catatan:"
  )

Ringkasan Koefisien Model dengan Ordinal sebagai Nominal

	Coefficients
Term	Estimate	Std. Error	Statistic	p-value
(Intercept)	-1.934213	0.579421	-3.3382	0.0008433
genderMale	-0.060080	0.269503	-0.2229	0.8235912
tech_savvinessBeginner	0.849712	0.353386	2.4045	0.0161953
tech_savvinessExpert	1.458630	0.474324	3.0752	0.0021038
tech_savvinessIntermediate	2.108440	0.343077	6.1457	0.0000000
income	0.037698	0.010164	3.7089	0.0002082
Catatan:
Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

Term	Estimate	Std. Error	z-value	p-value
(Intercept)	-1.93421	0.57942	-3.338	0.000843 ***
genderMale	-0.06008	0.26950	-0.223	0.823591
tech_savvinessBeginner	0.84971	0.35339	2.404	0.016195 *
tech_savvinessExpert	1.45863	0.47432	3.075	0.002104 **
tech_savvinessIntermediate	2.10844	0.34308	6.146	<0.000001 ***
income	0.03770	0.01016	3.709	0.000208 ***

Interpretasi :

1. Perlakuan Variabel sebagai Nominal (Dummy)

• Gender:
  Laki-laki (genderMale) memiliki peluang sukses lebih rendah dibanding perempuan, meskipun tidak signifikan secara statistik (*p-value = 0.7636*).

• Tech Savviness (Kategori):

  • Beginner: Signifikan meningkatkan peluang sukses (*Est. = 0.849, p = 0.016*).

  • Expert dan Intermediate: Efek lebih kuat dengan peningkatan peluang sukses yang sangat signifikan (p < 0.01).

• Income:
  Setiap kenaikan pendapatan meningkatkan peluang sukses secara signifikan (*Est. = 0.038, p < 0.001*).

library(knitr)
library(kableExtra)

# Buat data frame dari hasil summary
results <- data.frame(
  Term = c("(Intercept)", "genderMale", "tech_savviness_numeric", "income"),
  Estimate = c(0.146760, -0.075699, -0.256533, 0.029506),
  Std_Error = c(0.568490, 0.251707, 0.138778, 0.009237),
  z_value = c(0.258, -0.301, -1.849, 3.194),
  p_value = c(0.7963, 0.7636, 0.0645, 0.0014)
)

# Buat tabel
results %>%
  kable(format = "html", digits = 4,
        caption = "Ringkasan Koefisien Model Regresi Logistik",
        col.names = c("Term", "Estimate", "Std. Error", "z-value", "p-value")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  add_header_above(c(" ", "Coefficients" = 4)) %>%
  footnote(general = "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1",
           general_title = "Catatan:")

Ringkasan Koefisien Model Regresi Logistik

	Coefficients
Term	Estimate	Std. Error	z-value	p-value
(Intercept)	0.1468	0.5685	0.258	0.7963
genderMale	-0.0757	0.2517	-0.301	0.7636
tech_savviness_numeric	-0.2565	0.1388	-1.849	0.0645
income	0.0295	0.0092	3.194	0.0014
Catatan:
Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

Term	Estimate	Std. Error	z-value	p-value
(Intercept)	0.1468	0.5685	0.258	0.7963
genderMale	-0.0757	0.2517	-0.301	0.7636
tech_savviness_numeric	-0.2565	0.1388	-1.849	0.0645 .
income	0.0295	0.0092	3.194	0.0014 **

Interpretasi :

Perlakuan Variabel sebagai Numerik

• Tech Savviness (Skala Numerik):
  Setiap peningkatan satu unit skor tech savviness mengurangi peluang sukses (*Est. = -0.257*), tetapi efek ini hanya mendekati signifikansi (*p = 0.064*).

• Income:
  Konsisten signifikan dengan efek positif (*Est. = 0.030, p = 0.001*).

Pemilihan Model Regresi Logistik dan Evaluasi

Membangun Model Regresi Logistik : Pendekatan Confirmatory 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 Confrmatory dan Exploratory.

1. Confrmatory (Pendekatan Konfrmatori)

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 signifkan, 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 signifkan.

• Uji signifkansi 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 signifkan.

Explatory (Pendekatan Eksplatori)

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 signifkan.

• Backward Elimination: Mulai dari model penuh, variabel yang tidak signifkan 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 Confrmatory dan Exploratory bergantung pada tujuan penelitian. Jika ingin menguji hipotesis tertentu, gunakan pendekatan Confrmatory. Jika ingin menemukan model terbaik berdasarkan data, gunakan pendekatan Exploratory. Dalam praktiknya, kedua pendekatan ini sering digunakan secara komplementer: teori digunakan sebagai dasar, dan seleksi eksploratori dilakukan untuk menyempurnakan model.

Simulasi Data

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

## 
## Attaching package: 'MASS'

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

library(caret)

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

## Loading required package: lattice

## 
## Attaching package: 'caret'

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

library(pROC)

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

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

## 
## Attaching package: 'pROC'

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

library(DescTools)
set.seed(33)
n <- 233
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
lin_pred <- -0.5 + 1.2 * x1 - 0.8 * x2 + 0.5 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
df <- data.frame(y = as.factor(y), x1, x2, x3)
head(df)

Pemilihan Model

Model full

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

## 
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.9203     0.2604  -3.534  0.00041 ***
## x1            1.6384     0.2397   6.836 8.13e-12 ***
## x2           -0.7023     0.3543  -1.982  0.04746 *  
## x3            0.2300     0.1930   1.192  0.23326    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 286.50  on 232  degrees of freedom
## Residual deviance: 202.67  on 229  degrees of freedom
## AIC: 210.67
## 
## Number of Fisher Scoring iterations: 5

Metode Stepwise: Forward, Backward, dan Kedua Arah

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

Evaluasi Model: ROC dan AUC

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

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

## Setting direction: controls < cases

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

auc(roc_obj)

## Area under the curve: 0.8401

Pseudo R-Squared

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

##   CoxSnell Nagelkerke   McFadden 
##  0.2978826  0.4209776  0.2876123

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 149  31
##          1  13  40
##                                           
##                Accuracy : 0.8112          
##                  95% CI : (0.7549, 0.8593)
##     No Information Rate : 0.6953          
##     P-Value [Acc > NIR] : 4.241e-05       
##                                           
##                   Kappa : 0.5202          
##                                           
##  Mcnemar's Test P-Value : 0.01038         
##                                           
##             Sensitivity : 0.5634          
##             Specificity : 0.9198          
##          Pos Pred Value : 0.7547          
##          Neg Pred Value : 0.8278          
##              Prevalence : 0.3047          
##          Detection Rate : 0.1717          
##    Detection Prevalence : 0.2275          
##       Balanced Accuracy : 0.7416          
##                                           
##        'Positive' Class : 1               
## 

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

## Sensitivity Specificity 
##   0.5633803   0.9197531

Metode Perbandingan Model dalam Regresi Logistik

Tujuan

Dokumen ini menyajikan cara membandingkan model regresi logistik menggunakan ukuran :

- Deviance

- AIC (Akaike Information Criterion)

- Likelihood

- Ratio serta menjelaskan prinsip Parsimony dalam pemilihan model.

library(MASS)
library(broom)
library(DescTools)
set.seed(33)
n <- 333
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
lin_pred <- -1 + 1.2 * x1 - 0.6 * x2 + 0.8 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)
model1 <- glm(y ~ x1, data = data, family = binomial)
model2 <- glm(y ~ x1 + x2, data = data, family = binomial)
model3 <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
model_comp <- data.frame(
Model = c("Model 1", "Model 2", "Model 3"),
AIC = c(AIC(model1), AIC(model2), AIC(model3)),
Deviance = c(deviance(model1), deviance(model2), deviance(model3))
)
model_comp

Likelihood Ratio Test

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

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

Prinsip Parsimony

Model Selection Principles

Model yang kompleks sering memiliki AIC dan deviance yang lebih kecil. Namun:

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

Rumus dan Penjelasan

Rumus AIC \[ AIC = -2(\log L - k) = -2\log L + 2k \]

Penjelasan:
AIC adalah ukuran untuk menilai model berdasarkan kombinasi antara goodness-of-fit (melalui log-likelihood) dan kompleksitas (melalui jumlah parameter \(k\)). Semakin kecil AIC, semakin baik model tersebut secara keseluruhan karena AIC menghukum model yang terlalu kompleks meskipun memiliki likelihood tinggi.

Rumus Deviance

\[ D = -2\left[\log L(\text{model}) - \log L(\text{model saturasi})\right] \]

Penjelasan:
Deviance mengukur seberapa jauh model saat ini dibandingkan dengan model sempurna (model saturasi). Nilai deviance yang kecil menunjukkan model memberikan prediksi yang mendekati data aktual.

Rumus Likelihood-Ratio

\[ G^{2} = -2(\log L_{0} - \log L_{1}) \]

Penjelasan:
Statistik Likelihood Ratio digunakan untuk menguji apakah penambahan variabel dalam model secara signifikan meningkatkan kecocokan model. Jika \(G^{2}\) besar dan p-value kecil, maka model kompleks lebih baik dari model sederhana secara statistik.

Evaluasi Tabel Klasifikasi dan Akurasi Model

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

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 213  47
##          1  24  49
##                                           
##                Accuracy : 0.7868          
##                  95% CI : (0.7388, 0.8296)
##     No Information Rate : 0.7117          
##     P-Value [Acc > NIR] : 0.001174        
##                                           
##                   Kappa : 0.4405          
##                                           
##  Mcnemar's Test P-Value : 0.009030        
##                                           
##             Sensitivity : 0.5104          
##             Specificity : 0.8987          
##          Pos Pred Value : 0.6712          
##          Neg Pred Value : 0.8192          
##              Prevalence : 0.2883          
##          Detection Rate : 0.1471          
##    Detection Prevalence : 0.2192          
##       Balanced Accuracy : 0.7046          
##                                           
##        'Positive' Class : 1               
## 

Sensitivitas dan Spesifisitas

Sensitivitas

Kemampuan model mendeteksi kelas positif secara benar (True Positive Rate)

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

Spesifisitas

Kemampuan model mendeteksi kelas negatif secara benar (True Negative Rate)

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

---

Keterangan:
- TP (True Positive): Kasus positif yang terprediksi benar
- FN (False Negative): Kasus positif yang terprediksi salah sebagai negatif
- TN (True Negative): Kasus negatif yang terprediksi benar
- FP (False Positive): Kasus negatif yang terprediksi salah sebagai positif

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

## Sensitivity Specificity 
##   0.5104167   0.8987342

Kesimpulan Evaluasi Model

• Deviance: Nilai deviance yang kecil menunjukkan kecocokan model yang lebih baik dengan data observasi.
  Rumus:
  \[ D = -2(\log L_{\text{model}} - \log L_{\text{saturasi}}) \]

• AIC (Akaike Information Criterion): Nilai AIC yang rendah menunjukkan keseimbangan optimal antara kecocokan model dan kompleksitas (jumlah parameter).
  Rumus:
  \[ AIC = 2k - 2\log(L) \]

• Likelihood Ratio Test: Menguji signifikansi perbedaan antara model kompleks dan model sederhana.
  Rumus:
  \[ G^2 = -2(\log L_0 - \log L_1) \]

• Tabel Klasifikasi: Alat untuk mengevaluasi kinerja prediktif model dengan membandingkan hasil prediksi terhadap nilai aktual.
  Komponen utama:

  • True Positive (TP)
    
  • True Negative (TN)
    
  • False Positive (FP)
    
  • False Negative (FN)

• Prinsip Parsimoni (Occam’s Razor): Merekomendasikan pemilihan model yang lebih sederhana ketika performanya tidak berbeda signifikan dengan model yang lebih kompleks.

Detail ROC Penjelasan Kurva ROC (Receiver Operating Characteristic)

Kurva ROC (Receiver Operating Characteristic) adalah alat visual untuk mengevaluasi performa model klasifikasi biner yang menunjukkan trade-off antara:

• True Positive Rate (Sensitivity)
• False Positive Rate (1 - Specificity)

pada berbagai titik cut-off klasifikasi.

1. Definisi Komponen

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

Sumbu X:
\[ 1 - \text{Specificity} = \frac{FP}{FP + TN} \]

Karakteristik: - Garis diagonal menunjukkan performa acak (AUC = 0.5) - Kurva ideal mendekati pojok kiri atas (AUC ≈ 1)

2. Pengaruh Perubahan Cut-off

Saat cut-off diturunkan: - Model lebih mudah memprediksi positif - ➠ Sensitivitas meningkat (lebih banyak TP) - ➠ Spesifisitas menurun (lebih banyak FP)

Saat cut-off dinaikkan: - Model lebih ketat dalam memprediksi positif - ➠ Sensitivitas menurun (lebih banyak FN) - ➠ Spesifisitas meningkat (lebih banyak TN)

3. Bentuk Kurva Ideal

1. Awal vertikal tajam (mencapai sensitivitas tinggi cepat)
2. Kemudian horizontal sempurna (minim false positif)
3. AUC = 1 (luas area maksimal)

4. Interpretasi AUC

Nilai AUC	Klasifikasi Kinerja
0.90-1.00	Sangat Istimewa
0.80-0.89	Baik
0.70-0.79	Cukup
0.60-0.69	Marginal
0.50-0.59	Gagal

AUC dikenal juga sebagai concordance index, yaitu probabilitas bahwa model memberikan nilai skor probabilitas yang lebih tinggi untuk kasus positif daripada kasus negatif.

5. Kegunaan Kurva ROC

• Untuk membandingkan performa beberapa model klasifikasi
• Untuk memilih threshold (cut-off) optimal berdasarkan kebutuhan aplikasi (misalnya: lebih penting menghindari false negative atau false positive?)

6. Visualisasi dalam R

Kurva ROC dapat dibuat menggunakan package pROC :

library(pROC)
set.seed(33)
x1 <- rnorm(233)
x2 <- rbinom(233, 1, 0.5)
x3 <- rnorm(233)
lin_pred <- -1 + 1.5 * x1 - 0.7 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(233, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
pred <- predict(model, type = "response")
roc_obj <- roc(data$y, pred)

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

## Setting direction: controls < cases

plot(roc_obj)

auc(roc_obj)

## Area under the curve: 0.8794

7. Simulasi Pemilihan Threshold Optimal

Untuk memilih threshold terbaik, kita bisa mengevaluasi sensitivitas dan spesifsitas pada berbagai cut-of.

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

##    Threshold Sensitivity Specificity
## 1       0.10   0.9710145   0.5060976
## 2       0.15   0.8985507   0.6036585
## 3       0.20   0.8695652   0.6829268
## 4       0.25   0.8550725   0.7743902
## 5       0.30   0.8115942   0.8292683
## 6       0.35   0.7536232   0.8536585
## 7       0.40   0.6956522   0.8658537
## 8       0.45   0.6666667   0.9024390
## 9       0.50   0.5942029   0.9085366
## 10      0.55   0.5652174   0.9329268
## 11      0.60   0.4782609   0.9451220
## 12      0.65   0.4347826   0.9573171
## 13      0.70   0.3768116   0.9634146
## 14      0.75   0.3623188   0.9756098
## 15      0.80   0.3333333   0.9817073
## 16      0.85   0.2318841   0.9878049
## 17      0.90   0.1449275   1.0000000

Cut-of optimal bisa dipilih berdasarkan:

\(Maksimum dari Sensitivity + Specifcity\) Atau mempertimbangkan trade-of sesuai tujuan aplikasi (misalnya: jika False Negative harus dihindari, maka prioritaskan sensitivitas tinggi)

8. Catatan

ROC cocok saat proporsi kelas seimbang.

Untuk data dengan kelas tidak seimbang, precision-recall curve bisa lebih informatif.

Precision-Recall Curve (PR Curve)

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

1. Definisi

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

\[ \text{Precision} = \frac{TP}{TP + FP} \]

• Recall (Sensitivitas): Proporsi kasus positif yang berhasil diprediksi positif

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

2. Interpretasi

• PR Curve menunjukkan bagaimana presisi berubah saat recall meningkat.
• Idealnya, kita ingin nilai presisi dan recall keduanya tinggi, tetapi biasanya ada trade-off.
• Model dengan performa baik memiliki PR Curve yang melengkung ke pojok kanan atas.

3. Area Under PR Curve

• Luas kurva (AUPRC) mendekati 1 berarti model sangat baik.
• Baseline AUPRC = prevalensi kelas positif dalam data.

4. PR Curve vs ROC Curve

Aspek	ROC Curve	Precision-Recall Curve
Fokus	Semua kelas	Kelas positif saja
Kuat di	Data seimbang	Data tidak seimbang
Sumbu Y	Sensitivitas (Recall)	Precision
Sumbu X	1 - Spesifisitas	Recall

5. Visualisasi PR Curve di R

library(PRROC)

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

## Loading required package: rlang

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

set.seed(33)
x1 <- rnorm(233)
x2 <- rbinom(233, 1, 0.5)
x3 <- rnorm(233)
lin_pred <- -1 + 1.5 * x1 - 0.7 * x2 + 0.6 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(233, 1, p)
data <- data.frame(y = y, x1, x2, x3)
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
prob <- predict(model, type = "response")
pr <- pr.curve(scores.class0 = prob[data$y == 1],
scores.class1 = prob[data$y == 0],
curve = TRUE)
plot(pr)

6. Catatan

• PR Curve sangat informatif untuk aplikasi seperti deteksi penipuan atau diagnosis penyakit langka.

• Gunakan PR Curve saat:

  • Kelas positif jauh lebih jarang daripada kelas negatif

  • Tujuan aplikasi lebih mementingkan presisi terhadap kelas minoritas

Pseudo R-squared pada Regresi Logistik

Tujuan

Dokumen ini menjelaskan dan menghitung pseudo R-Squared dalam regresi logistik :

• \(R^2_{\text{CoxandSnell}}\)
• \(R^2_{\text{McFadden}}\)

Simulasi Data

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

Model logistik dan Null model

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

Likelihood dan Rumus

\[ R^2_{\text{Cox and Snell}} = 1 - \left(\frac{L_0}{L_M}\right)^{2/n} \]

\[ R^2_{\text{McFadden}} = 1 - \frac{\log L_M}{\log L_0} \]

Dengan:

• \(L_0\): likelihood model null (tanpa prediktor)
• \(L_M\): likelihood model penuh

Perhitungan Manual R-squared

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

Perhitungan Otomatis dengan Package Tambahan

Menggunakan pacl

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

## Loading required package: pscl

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

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

pR2(model)

## fitting null model for pseudo-r2

##          llh      llhNull           G2     McFadden         r2ML         r2CU 
## -144.0184091 -200.0037693  111.9707205    0.2799215    0.2855544    0.4084167

Menggunakan desctools

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

##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.2799215       0.2599219       0.2855544       0.4084167       0.2516362 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.4611193       0.3154733       0.4560219       0.3184918     296.0368181 
##             BIC          logLik         logLik0              G2 
##     311.2693881    -144.0184091    -200.0037693     111.9707205

Interpretasi

• Nilai \(R^2\) mendekati 1 berarti model memiliki kekuatan prediktif yang baik.

• McFadden \(R^2 > 0.2\) sering dianggap sebagai model dengan kecocokan yang baik.

• Cox & Snell \(R^2\) lebih konservatif dan tidak pernah mencapai 1 penuh.

Gunakan beberapa pendekatan ini sebagai pelengkap untuk menilai performa model logistik secara lebih menyeluruh.

Apa itu Distribusi Multinomial

Distribusi multinomial adalah perluasan dari distribusi binomial untuk lebih dari dua kategori. Jika \(X_1, X_2, \ldots, X_k\) menyatakan banyaknya kejadian dalam masing-masing dari \(k\) 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\) dan \(\sum_{i=1}^k p_i = 1\).

Studi Kasus

Sebuah survei dilakukan terhadap 20 orang yang diminta memilih salah satu dari tiga jenis sepatu favorit: Boots (B), Sneakers (S), dan Heels (H).

Hasil survei:

• Boots: 6 orang
• Sneakers: 8 orang
• Heels: 6 orang

Probabilitas teoritik preferensi:

• \(p_B = 0.6\),
• \(p_S = 0.8\),
• \(p_H = 0.6\)

Pertanyaannya: Berapa peluang bahwa dalam 20 orang akan ada 6 yang memilih boots, 8 memilih sneakers, dan 6 memilih heels?

Rumus Distribusi Multinomial Distribusi peluang multinomial:

\[ 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:

• \(n = 10\), \(x_1 = 4\), \(x_2 = 3\), \(x_3 = 3\)
• \(p_1 = 0.3\), \(p_2 = 0.4\), \(p_3 = 0.3\)

Perhitungan Manual di R

n <- 20
x <- c(6, 8, 6)
p <- c(0.3, 0.4, 0.3)
# Hitung komponen-koefisien
faktorial_total <- factorial(n)
faktorial_x <- prod(factorial(x))
koefisien <- faktorial_total / faktorial_x
# Hitung peluang
peluang <- koefisien * prod(p^x)
peluang

## [1] 0.0405391

Probabilitas bahwa 6 orang memilih boots, 8 sneakers, dan 6 heels (dengan proporsi preferensi 0.3, 0.4, dan 0.3) adalah 0.0405391. Distribusi multinomial digunakan untuk menghitung peluang dalam percobaan dengan beberapa kategori hasil. Rumus dasarnya merupakan generalisasi dari binomial untuk lebih dari dua kategori.

Multinomial Logistik Regression

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

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

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

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

Baseline Category Logit Model

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

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

dengan:

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

Maka, terdapat sebanyak \((c-1)\) fungsi logit.

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

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

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

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

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

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

Terdapat dua model logit, satu untuk perbandingan kategori 1 dengan 3, dan satu lagi untuk kategori 2 dengan 3.

Relasi Antar Kategori Jika ingin menghitung logit antara kategori 1 dan 2:

\[ \log \left(\frac{\pi_1}{\pi_2}\right) = \log \left(\pi_1 / \pi_3\right) - \log \left(\pi_2 / \pi_3\right) \]

\[ = (\alpha_1 + \beta_1 x) - (\alpha_2 + \beta_2 x) = (\alpha_1 - \alpha_2) + (\beta_1 - \beta_2)x \]

Model baseline-category logit:

• Digunakan untuk respon dengan kategori \(> 2\)
• Menghasilkan \((c-1)\) fungsi logit terhadap satu baseline
• Logit antara kategori selain baseline dapat dihitung dari selisih dua logit terhadap baseline

Implementasi di R menggunakan fungsi multinom() dari package nnet, dan kategori baseline bisa diatur dengan relevel().

Estimasi Parameter

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

Log-likelihood:

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

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

Contoh Kasus

Sebuah penyedia layanan streaming video ingin memahami faktor-faktor yang memengaruhi retensi pelanggan mereka, khususnya jenis paket langganan yang dipilih pelanggan setelah masa percobaan gratis berakhir. Perusahaan mencurigai bahwa durasi penggunaan selama masa percobaan, genre film/serial favorit, dan interaksi dengan fitur platform dapat memengaruhi keputusan langganan.

Perusahaan melakukan survei terhadap 300 pengguna yang baru saja menyelesaikan masa percobaan gratis dan mengumpulkan data berikut:

• Keputusan Langganan: Jenis paket yang dipilih (Tidak Berlangganan, Paket Dasar, Paket Premium).
• Durasi Penggunaan Percobaan: Total jam menonton selama masa percobaan gratis.
• Genre Favorit: Genre yang paling banyak ditonton (misalnya, Drama, Komedi, Aksi).
• Interaksi Fitur: Indeks aktivitas pengguna dengan fitur tambahan (misalnya, membuat daftar putar, menggunakan fitur skip intro, mengunduh konten).

Tujuannya adalah:

Mengetahui bagaimana durasi penggunaan masa percobaan, genre favorit, dan interaksi fitur memengaruhi pilihan paket langganan pelanggan.

Simulasi Data

set.seed(123) # Untuk reproduktibilitas hasil
n <- 300 # Jumlah pelanggan yang disurvei, sesuai contoh kasus

# Variabel Prediktor
Usia_Pelanggan <- round(rnorm(n, mean = 30, sd = 8)) # Usia Pelanggan (misalnya 18-60 tahun)
Durasi_Penggunaan_Percobaan <- round(pmax(rnorm(n, mean = 10, sd = 4), 0)) # Durasi penggunaan dalam jam, minimal 0

# --- Perubahan di sini: Hanya sertakan genre yang diinginkan ---
Genre_Favorit <- sample(c("Drama", "Komedi", "Aksi"), n, replace = TRUE) # Hanya Drama, Komedi, Aksi
# --- Akhir perubahan ---

Interaksi_Fitur <- round(rnorm(n, mean = 5, sd = 2)) # Indeks interaksi fitur (misalnya 0-10)
Interaksi_Fitur <- pmax(Interaksi_Fitur, 0) # Pastikan tidak ada nilai negatif

# Simulasikan Keputusan_Langganan berdasarkan probabilitas yang bervariasi
Keputusan_Langganan <- sapply(1:n, function(i) {
  prob_tidak_berlangganan <- 0.4 # Baseline probabilitas
  prob_paket_dasar <- 0.3
  prob_paket_premium <- 0.3

  # Sesuaikan probabilitas berdasarkan beberapa prediktor (contoh sederhana)
  if (Durasi_Penggunaan_Percobaan[i] > 10) {
    prob_paket_dasar <- prob_paket_dasar + 0.1
    prob_paket_premium <- prob_paket_premium + 0.1
    prob_tidak_berlangganan <- prob_tidak_berlangganan - 0.2
  }
  # Sesuaikan logika ini jika Anda ingin probabilitas berbeda untuk genre spesifik
  if (Genre_Favorit[i] == "Aksi") {
    prob_paket_premium <- prob_paket_premium + 0.08
    prob_paket_dasar <- prob_paket_dasar - 0.03
    prob_tidak_berlangganan <- prob_tidak_berlangganan - 0.05
  } else if (Genre_Favorit[i] == "Drama") {
    prob_paket_dasar <- prob_paket_dasar + 0.05
    prob_tidak_berlangganan <- prob_tidak_berlangganan - 0.05
  } # Komedi tidak ada penyesuaian khusus di sini, akan menggunakan base probabilitas

  if (Interaksi_Fitur[i] > 7) {
    prob_paket_premium <- prob_paket_premium + 0.1
    prob_paket_dasar <- prob_paket_dasar + 0.05
    prob_tidak_berlangganan <- prob_tidak_berlangganan - 0.15
  }

  # Normalisasi probabilitas agar totalnya 1
  total_prob <- prob_tidak_berlangganan + prob_paket_dasar + prob_paket_premium
  # Pastikan tidak ada probabilitas negatif yang bisa muncul dari penyesuaian di atas
  probs_adjusted <- pmax(c(prob_tidak_berlangganan, prob_paket_dasar, prob_paket_premium), 0)
  probs <- probs_adjusted / sum(probs_adjusted)

  sample(c("Tidak Berlangganan", "Paket Dasar", "Paket Premium"), size = 1, prob = probs)
})

# Buat dataframe
df_ecommerce <- data.frame(
  Keputusan_Langganan = factor(Keputusan_Langganan),
  Usia_Pelanggan = Usia_Pelanggan,
  Durasi_Penggunaan_Percobaan = Durasi_Penggunaan_Percobaan,
  Genre_Favorit = factor(Genre_Favorit), # Penting untuk tetap sebagai faktor
  Interaksi_Fitur = Interaksi_Fitur
)

# Set 'Tidak Berlangganan' sebagai kategori referensi (baseline)
df_ecommerce$Keputusan_Langganan <- relevel(df_ecommerce$Keputusan_Langganan, ref = "Tidak Berlangganan")

# Verifikasi level faktor Genre_Favorit
print(levels(df_ecommerce$Genre_Favorit))

## [1] "Aksi"   "Drama"  "Komedi"

# Tampilkan beberapa baris pertama dari dataframe
head(df_ecommerce)

Estimasi Model

# Kemudian, jalankan model multinom
library(nnet)
model_mnlogit_ecommerce <- multinom(Keputusan_Langganan ~ Usia_Pelanggan + Durasi_Penggunaan_Percobaan + Genre_Favorit + Interaksi_Fitur, data = df_ecommerce)

## # weights:  21 (12 variable)
## initial  value 329.583687 
## iter  10 value 308.639676
## final  value 307.598261 
## converged

summary(model_mnlogit_ecommerce)

## Call:
## multinom(formula = Keputusan_Langganan ~ Usia_Pelanggan + Durasi_Penggunaan_Percobaan + 
##     Genre_Favorit + Interaksi_Fitur, data = df_ecommerce)
## 
## Coefficients:
##               (Intercept) Usia_Pelanggan Durasi_Penggunaan_Percobaan
## Paket Dasar     -1.029037   -0.007663560                   0.1567823
## Paket Premium   -1.222637    0.007327071                   0.1364632
##               Genre_FavoritDrama Genre_FavoritKomedi Interaksi_Fitur
## Paket Dasar            0.3492730          -0.1281806     0.004600098
## Paket Premium         -0.4059777          -0.6631106     0.120957818
## 
## Std. Errors:
##               (Intercept) Usia_Pelanggan Durasi_Penggunaan_Percobaan
## Paket Dasar     0.8838777     0.02088255                  0.04214324
## Paket Premium   0.8585725     0.02016313                  0.04048041
##               Genre_FavoritDrama Genre_FavoritKomedi Interaksi_Fitur
## Paket Dasar            0.4045929           0.3979505      0.07710452
## Paket Premium          0.3873188           0.3728600      0.07427517
## 
## Residual Deviance: 615.1965 
## AIC: 639.1965

Nilai P-Value dan Interpretasi

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

##               (Intercept) Usia_Pelanggan Durasi_Penggunaan_Percobaan
## Paket Dasar        0.2443         0.7136                       2e-04
## Paket Premium      0.1544         0.7163                       7e-04
##               Genre_FavoritDrama Genre_FavoritKomedi Interaksi_Fitur
## Paket Dasar               0.3880              0.7474          0.9524
## Paket Premium             0.2946              0.0753          0.1034

Interpretasi :

• Koefsien untuk kategori “Desktop” dan “Tablet” dibandingkan dengan baseline “Laptop”

• Nilai p-value kecil (<0.05) menunjukkan variabel tersebut signifkan memengaruhi preferensi perangkat

Prediksi dan Validasi

# --- Prediction Code (Corrected) ---
# Add predicted class to the dataframe
df_ecommerce$Predicted_Keputusan_Langganan <- predict(model_mnlogit_ecommerce, newdata = df_ecommerce, type = "class")

# Add predicted probabilities for each class to the dataframe
# This will create multiple columns, one for each probability (e.g., Tidak.Berlangganan, Paket.Dasar, Paket.Premium)
predicted_probs <- predict(model_mnlogit_ecommerce, newdata = df_ecommerce, type = "probs")
df_ecommerce <- cbind(df_ecommerce, predicted_probs) # Combine with original dataframe

# You can also view the table of actual vs. predicted classes
cat("\nConfusion Matrix (Actual vs. Predicted):\n")

## 
## Confusion Matrix (Actual vs. Predicted):

print(table(Actual = df_ecommerce$Keputusan_Langganan, Predicted = df_ecommerce$Predicted_Keputusan_Langganan))

##                     Predicted
## Actual               Tidak Berlangganan Paket Dasar Paket Premium
##   Tidak Berlangganan                 20          15            38
##   Paket Dasar                        13          38            51
##   Paket Premium                      15          30            80

Kesimpulan

Model regresi logistik multinomial berhasil digunakan untuk :

• Menganalisis hubungan antara atribut karyawan dan preferensi perangkat kerja

• Mengetahui faktor signifkan yang memengaruhi pilihan

• Memungkinkan prediksi jenis perangkat yang dipilih oleh karyawan baru berdasarkan karakteristiknya

Contoh Kasus 2 di R

Kita gunakan dataset iris untuk memodelkan jenis spesies bungan berdasarkan lebar dan panjang kelopak.

library(dplyr)
data(iris)
iris <- iris %>% mutate(Species = relevel(Species, ref = "setosa"))
model <- multinom(Species ~ Petal.Length + Petal.Width, data = iris)

## # weights:  12 (6 variable)
## initial  value 164.791843 
## iter  10 value 12.657828
## iter  20 value 10.374056
## iter  30 value 10.330881
## iter  40 value 10.306926
## iter  50 value 10.300057
## iter  60 value 10.296452
## iter  70 value 10.294046
## iter  80 value 10.292029
## iter  90 value 10.291154
## iter 100 value 10.289505
## final  value 10.289505 
## stopped after 100 iterations

summary(model)

## Call:
## multinom(formula = Species ~ Petal.Length + Petal.Width, data = iris)
## 
## Coefficients:
##            (Intercept) Petal.Length Petal.Width
## versicolor   -22.79944      6.92122    7.878496
## virginica    -67.82521     12.64721   18.261016
## 
## Std. Errors:
##            (Intercept) Petal.Length Petal.Width
## versicolor     44.3859     37.58715    81.00888
## virginica      46.3939     37.65702    81.09482
## 
## Residual Deviance: 20.57901 
## AIC: 32.57901

Interpretasi Koefisien

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

##            (Intercept) Petal.Length Petal.Width
## versicolor      0.6075       0.8539      0.9225
## virginica       0.1438       0.7370      0.8218

Nilai p-value menunjukkan apakah variabel prediktor berpengaruh signifkan terhadap log odds dibanding baseline category.

Prediksi dan Visualisasi

iris$predicted <- predict(model, newdata = iris)
table(Predicted = iris$predicted, Actual = iris$Species)

##             Actual
## Predicted    setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         47         3
##   virginica       0          3        47

ggplot(iris, aes(x = Petal.Length, y = Petal.Width, color = predicted)) +
geom_point(size = 2) +
labs(title = "Multinomial Logistic Regression Predictions",
x = "Petal Length", y = "Petal Width") +
theme_minimal()

Kesimpulan Multinomial logistic regression efektif digunakan untuk klasifkasi kategori nominal. Model ini memberikan estimasi log-odds terhadap baseline dan prediksi kategori baru berdasarkan kombinasi prediktor.

Regresi Logistik Ordinal

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

Model ini berbeda dengan: * Regresi logistik biner: hanya 2 kategori * Regresi logistik multinomial: kategori \(> 2\) tetapi tidak berurutan

Konsep Cumulative Logit Model

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

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

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

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

Interpretasi Koefisien

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

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

\[ OR = e^\beta \]

Contoh Data: Kepuasan Pelanggan

Misal kita memiliki data fiktif tingkat kepuasan pelanggan (1: Tidak Puas, 2: Cukup, 3: Puas) terhadap kecepatan layanan:

set.seed(123)
n <- 200
speed <- round(runif(n, 1, 10))
satisfaction <- cut(5 + 0.5*speed + rnorm(n),
breaks = c(-Inf, 5.5, 7.5, Inf),
labels = c("Tidak Puas", "Cukup", "Puas"),
ordered_result = TRUE)
df <- data.frame(satisfaction, speed)
head(df)

Estimasi Model Ordinal

model_ord <- polr(satisfaction ~ speed, data = df, Hess = TRUE)
summary(model_ord)

## Call:
## polr(formula = satisfaction ~ speed, data = df, Hess = TRUE)
## 
## Coefficients:
##        Value Std. Error t value
## speed 0.9096     0.1094   8.315
## 
## Intercepts:
##                  Value  Std. Error t value
## Tidak Puas|Cukup 1.3015 0.4377     2.9738 
## Cukup|Puas       4.4734 0.5718     7.8232 
## 
## Residual Deviance: 237.2312 
## AIC: 243.2312

Nilai P-Value

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

##                      Value Std. Error  t value
## speed            0.9095585  0.1093925 8.314630
## Tidak Puas|Cukup 1.3015075  0.4376597 2.973789
## Cukup|Puas       4.4733938  0.5718127 7.823180

p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = round(p, 4)))

##                      Value Std. Error  t value p value
## speed            0.9095585  0.1093925 8.314630  0.0000
## Tidak Puas|Cukup 1.3015075  0.4376597 2.973789  0.0029
## Cukup|Puas       4.4733938  0.5718127 7.823180  0.0000

Prediksi Probabilitas

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

##    Tidak Puas     Cukup      Puas
## 1 0.037460604 0.4439482 0.5185912
## 2 0.015430723 0.2566765 0.7278928
## 3 0.006271788 0.1245723 0.8691559
## 4 0.002535158 0.0546231 0.9428417
## 5 0.001022461 0.0228089 0.9761686

Goodness of Fit dan Proportional Odds

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

Alternatif Model Ordinal

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

Model tersebut dapat digunakan saat asumsi proportional odds tidak terpenuhi.

Kesimpulan

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

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

Asumsi Paralelisme dalam Regresi Logistik Ordinal

Model regresi logistik ordinal yang paling umum digunakan adalah Cumulative Logit Model dengan asumsi Proportional Odds. Asumsi ini dikenal juga sebagai asumsi paralelisme (parallel lines assumption). Definisi Asumsi Paralelisme menyatakan bahwa koefisien regresi (\(\beta\)) sama untuk setiap kategori kumulatif dari variabel respon.

Bentuk umum model:

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

Untuk \(j = 1, \ldots, c-1\).

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

Visualisasi Dalam asumsi paralelisme, kurva logit kumulatif dari tiap kategori terhadap prediktor akan memiliki kemiringan yang sama (paralel), hanya berbeda pada posisi (intercept).

Konsekuensi Pelanggaran Asumsi Jika asumsi ini tidak terpenuhi:

• Efek prediktor berbeda untuk setiap batas kategori.
• Model cumulative logit tidak valid.
• Perlu pengunaan model alternatif:
  • Generalized Ordinal Logistic Regression
  • Partial Proportional Odds Model

Pengujian Asumsi Paralelisme Untuk memeriksa validitas asumsi, dapat digunakan:

• Likelihood Ratio Test antara model proporsional dan non-proporsional
• Brant Test (paket brant di R)

Contoh (dengan package brant):

# 1. Instal dan muat package yang dibutuhkan
# install.packages("MASS")
# install.packages("brant")
# install.packages("lmtest")
# install.packages("VGAM") # Diperlukan untuk model non-proporsional (vglm)

library(MASS)    # Untuk fungsi polr()
library(brant)   # Untuk fungsi brant()

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

library(lmtest)  # Untuk fungsi lrtest()    

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

## Loading required package: zoo

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

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

# 2. Simulasi Data (sesuai contoh kasus perangkat kerja)
set.seed(123)
n <- 150
# Asumsikan Device memiliki urutan, misalnya Laptop < Desktop < Tablet (atau sebaliknya)
# Ini penting untuk regresi logistik ordinal.
# Mari kita buat asumsi bahwa ada "preferensi" atau "urutan" yang tersirat.
# Contoh: Laptop (paling sederhana/mobile), Desktop (lebih kuat), Tablet (di tengah-tengah/spesifik)
# Untuk tujuan simulasi, kita bisa mengurutkan level secara manual.
Device_ordered <- sample(c("Laptop", "Tablet", "Desktop"), n, replace = TRUE,
                         prob = c(0.4, 0.3, 0.3)) # Sesuaikan probabilitas jika ingin bias

Department <- sample(c("IT", "Marketing", "HR"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 35, sd = 7))
Experience <- round(pmax(rnorm(n, mean = 7, sd = 3), 0))

# Simulasikan Device dengan efek prediktor agar ada hubungan
# Logika sederhana untuk menghasilkan data ordinal dengan prediktor
sim_Device <- function(age, dept, exp) {
  # Base probabilitas untuk setiap kategori
  p_laptop <- 0.3
  p_tablet <- 0.3
  p_desktop <- 0.4

  # Pengaruh Usia: cenderung ke Desktop jika lebih tua
  if (age > 40) {
    p_desktop <- p_desktop + 0.1
    p_laptop <- p_laptop - 0.1
  }
  # Pengaruh Departemen: IT cenderung ke Desktop/Laptop, Marketing ke Tablet/Laptop
  if (dept == "IT") {
    p_desktop <- p_desktop + 0.1
    p_tablet <- p_tablet - 0.1
  } else if (dept == "Marketing") {
    p_tablet <- p_tablet + 0.1
    p_desktop <- p_desktop - 0.1
  }
  # Pengaruh Pengalaman: lebih berpengalaman cenderung ke perangkat yang disesuaikan
  if (exp > 10) {
    p_desktop <- p_desktop + 0.05
    p_tablet <- p_tablet - 0.05
  }

  # Normalisasi probabilitas dan pastikan non-negatif
  probs <- pmax(c(p_laptop, p_tablet, p_desktop), 0)
  probs <- probs / sum(probs)

  sample(c("Laptop", "Tablet", "Desktop"), 1, prob = probs)
}

# Terapkan simulasi ke setiap baris data
Device_simulated <- mapply(sim_Device, Age, Department, Experience)

# Buat dataframe
df_perangkat_kerja <- data.frame(
  Device = factor(Device_simulated, levels = c("Laptop", "Tablet", "Desktop"), ordered = TRUE), # Penting: ordered = TRUE
  Age = Age,
  Department = factor(Department),
  Experience = Experience
)

cat("Head of df_perangkat_kerja:\n")

## Head of df_perangkat_kerja:

print(head(df_perangkat_kerja))

##    Device Age Department Experience
## 1 Desktop  35         HR          3
## 2  Tablet  50  Marketing          9
## 3  Tablet  30         HR          5
## 4 Desktop  27         IT          9
## 5  Tablet  35         HR          7
## 6 Desktop  37         HR          9

# 3. Fitting Model Regresi Logistik Ordinal (Cumulative Logit Model)
# Menggunakan polr() dari package MASS
cat("\nFitting Ordinal Logistic Regression Model (polr):\n")

## 
## Fitting Ordinal Logistic Regression Model (polr):

model_polr <- polr(Device ~ Age + Department + Experience, data = df_perangkat_kerja, Hess = TRUE)
print(summary(model_polr))

## Call:
## polr(formula = Device ~ Age + Department + Experience, data = df_perangkat_kerja, 
##     Hess = TRUE)
## 
## Coefficients:
##                        Value Std. Error t value
## Age                  0.02349    0.02375  0.9889
## DepartmentIT         0.55006    0.38771  1.4187
## DepartmentMarketing -0.43401    0.38425 -1.1295
## Experience          -0.05407    0.05240 -1.0320
## 
## Intercepts:
##                Value   Std. Error t value
## Laptop|Tablet  -0.5932  0.9794    -0.6057
## Tablet|Desktop  0.6435  0.9782     0.6579
## 
## Residual Deviance: 309.1045 
## AIC: 321.1045

# 4. Pengujian Asumsi Paralelisme

# Uji 1: Brant Test
cat("\nBrant Test for Parallelism Assumption:\n")

## 
## Brant Test for Parallelism Assumption:

brant_test_result <- brant(model_polr)

## ---------------------------------------------------- 
## Test for     X2  df  probability 
## ---------------------------------------------------- 
## Omnibus          16.55   4   0
## Age          0.56    1   0.45
## DepartmentIT     8.62    1   0
## DepartmentMarketing  0   1   0.99
## Experience       0.25    1   0.62
## ---------------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds

print(brant_test_result)

##                               X2 df probability
## Omnibus             1.655114e+01  4 0.002362142
## Age                 5.610983e-01  1 0.453818073
## DepartmentIT        8.615380e+00  1 0.003333362
## DepartmentMarketing 2.559945e-04  1 0.987234528
## Experience          2.451264e-01  1 0.620527853

# Interpretasi: Jika p-value < 0.05, maka asumsi paralelisme tidak terpenuhi.

Jika p-value < 0.05 \(\rightarrow\) asumsi tidak terpenuhi. Kesimpulan

• Asumsi paralelisme penting untuk validitas model cumulative logit.
• Menyederhanakan interpretasi karena efek prediktor konstan.
• Jika tidak terpenuhi, gunakan model ordinal alternatif.

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.

Pendekatan statistik yang umum digunakan dalam analisis data kategorik meliputi:

1. Tabel Kontingensi
   Digunakan untuk eksplorasi hubungan antar dua atau lebih variabel kategorik dengan cara menyajikan frekuensi gabungan. Tabel ini bisa dilengkapi dengan uji asosiasi seperti chi-square untuk mengetahui hubungan antar variabel.

2. Model Log-linear
   Digunakan untuk menganalisis struktur asosiasi dari tabel kontingensi tanpa variabel dependen. Model ini cocok untuk memahami interaksi antar variabel kategorik dan dapat diuji menggunakan likelihood ratio test.

3. Model Regresi Logistik
   Cocok digunakan saat ada variabel kategorik sebagai variabel dependen (contoh: sakit/tidak). Model ini memodelkan peluang terjadinya suatu kejadian (log odds) berdasarkan variabel independen dan dapat diperluas menjadi regresi logistik multinomial atau ordinal.

Perbandingan

• Tabel kontingensi: bersifat deskriptif dan eksploratif.
• Model log-linear: fokus pada struktur asosiasi variabel.
• Regresi logistik: fokus pada prediksi terhadap kategori outcome.

Kesimpulan

Pemilihan pendekatan tergantung pada tujuan analisis: apakah untuk deskripsi, eksplorasi struktur, atau prediksi. Kombinasi pendekatan sering digunakan untuk pemahaman yang lebih komprehensif terhadap data kategorik.

Tabel Kontingensi

Tabel kontingensi menyajikan jumlah frekuensi dari kombinasi kategori antar variabel.

Contoh:

table_data <- matrix(c(30, 20, 50, 70), nrow=2, 
       dimnames = list(Obat = c("Timolol", "Placebo"),
                       Serangan = c("Ya", "Tidak")))
table_data

##          Serangan
## Obat      Ya Tidak
##   Timolol 30    50
##   Placebo 20    70

Tabel kontingensi bersifat deskriptif dan tidak melibatkan pemodelan probabilitas.

Model Loglinear

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

\[ \log(\mu_{ij}) = \lambda + \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(~ Obat * Serangan, data = table_data)

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

Model Regresi Logistik

Model regresi logistik biner:

\[ \log\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(
  Serangan = c(1, 0, 1, 0),
  Obat = factor(c("Timolol", "Timolol", "Placebo", "Placebo")),
  Frek = c(30, 20, 50, 70)
)
model_logit <- glm(Serangan ~ Obat, weights = Frek, family = binomial, data = data_glm)
summary(model_logit)

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

Aspek	Tabel Kontingensi	Model Loglinear	Regresi Logistik
Tujuan	Deskripsi frekuensi	Deteksi asosiasi	Prediksi probabilitas
Variabel Dependen	Tidak ada	Tidak ada (simetris)	Ada (eksplisit)
Distribusi	Tidak diasumsikan	Poisson (frekuensi sel)	Binomial (probabilitas)
Bentuk Model	Tidak ada	GLM: log(μ) ~ efek	GLM: logit(p) ~ prediktor
Cocok untuk	Eksplorasi awal	Tabel > 2 variabel	Studi prediktif

Tabel Kontigensi dan Model Loglinier

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

# Contoh tabel 2x2
matrix(c(30, 20, 50, 70), nrow=2,
       dimnames = list(Obat = c("Timolol", "Placebo"),
                       Serangan = c("Ya", "Tidak")))

##          Serangan
## Obat      Ya Tidak
##   Timolol 30    50
##   Placebo 20    70

Model log-linier untuk tabel I x J dapat dituliskan: \[ \log(\mu_{ij}) = \lambda + \lambda_i^{T} + \lambda_j^{R} + \lambda_{ij}^{TR} \]

Model Saturated

Model saturated atau model penuh menyertakan seluruh efek utama dan interaksi:

• Cocok sempurna terhadap data
• Tidak mengasumsikan independensi antar variabel Contoh formulasi untuk tabel 2x2:

# Data
library(MASS)
data <- matrix(c(35, 65, 45, 55), nrow=2, byrow=TRUE)
dimnames(data) <- list(Obat = c("Timolol", "Placebo"), Serangan = c("Ya", "Tidak"))
ftable(data)

##         Serangan Ya Tidak
## Obat                     
## Timolol          35    65
## Placebo          45    55

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

model_saturated <- loglm(~ Obat * Serangan, data = data)
summary(model_saturated)

## Formula:
## ~Obat * Serangan
## attr(,"variables")
## list(Obat, Serangan)
## attr(,"factors")
##          Obat Serangan Obat:Serangan
## Obat        1        0             1
## Serangan    0        1             1
## attr(,"term.labels")
## [1] "Obat"          "Serangan"      "Obat:Serangan"
## 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}) = \lambda + \lambda_i^{T} + \lambda_j^{R} \]

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

model_indep <- loglm(~ Obat + Serangan, data = data)
summary(model_indep)

## Formula:
## ~Obat + Serangan
## attr(,"variables")
## list(Obat, Serangan)
## attr(,"factors")
##          Obat Serangan
## Obat        1        0
## Serangan    0        1
## attr(,"term.labels")
## [1] "Obat"     "Serangan"
## 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 2.087576  1 0.1485015
## Pearson          2.083333  1 0.1489147

Odds Ratio dan Interpretasi

Odds ratio untuk tabel 2x2:

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

Interpretasi nilai OR:

• OR = 1: Tidak ada asosiasi
• OR > 1: Asosiasi positif
• OR < 1: Asosiasi negatif

Estimasi Parameter

Dalam model saturated:

• Estimasi dilakukan dengan pembatasan seperti sum-to-zero
• Estimasi parameter dilakukan dengan iterative proportional fitting (IPF)

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

## [1] -0.4183685

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:
##  ~Obat + Serangan 
## Model 2:
##  ~Obat * Serangan 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   2.087576  1                                    
## Model 2   0.000000  0   2.087576         1         0.1485
## Saturated 0.000000  0   0.000000         0         1.0000

Studi Kasus: Kepercayaan terhadap Surga

Studi dari Agresti (2019) membahas hubungan antara kebahagiaan dan kepercayaan terhadap kehidupan akhirat.

data_survey <- matrix(c(32,190,
                        113,611,
                        51,326),
                      nrow = 3, byrow = TRUE,
                      dimnames = list(Kebahagiaan = c("Tidak", "Cukup", "Sangat"),
                                      Surga = c("Tidak Percaya", "Percaya")))
ftable(data_survey)

##             Surga Tidak Percaya Percaya
## Kebahagiaan                            
## Tidak                        32     190
## Cukup                       113     611
## Sangat                       51     326

loglm(~ Kebahagiaan + Surga, data = data_survey)

## Call:
## loglm(formula = ~Kebahagiaan + Surga, data = data_survey)
## 
## Statistics:
##                        X^2 df  P(> X^2)
## Likelihood Ratio 0.8911136  2 0.6404675
## Pearson          0.8836760  2 0.6428538

Rangkuman 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 (μij ) dapat dinyatakan sebagai penjumlahan efek variabel dan (bila perlu) interaksinya. Untuk tabel 2x2: \[ \log(\mu_{ij}) = \lambda + \lambda_i^{A} + \lambda_j^{B} + \lambda_{ij}^{AB} \]

Perbedaan Utama Antara Model Log Linear dan Model Regresi Logistik - Model log-linear

digunakan untuk memodelkan frekuensi (count) pada tabel kontingensi dan menguji asosiasi antar variabel kategorik, tanpa menganggap ada variabel respon dan prediktor.

• Model regresi logistik digunakan untuk memodelkan probabilitas kejadian suatu outcome (biner) berdasarkan satu atau lebih prediktor (bisa kategorik maupun kontinu).

Estimasi Parameter Log Linear 2 Arah

Sistem Persamaan Model Log-Linear

Sistem Persamaan Model Log-Linear

\[ \log(\mu_{11}) = \lambda + \lambda^A_1 + \lambda^B_1 + \lambda^{AB}_{11} \]

\[ \log(\mu_{12}) = \lambda + \lambda^A_1 + \lambda^B_2 + \lambda^{AB}_{12} \]

\[ \log(\mu_{21}) = \lambda + \lambda^A_2 + \lambda^B_1 + \lambda^{AB}_{21} \]

\[ \log(\mu_{22}) = \lambda + \lambda^A_2 + \lambda^B_2 + \lambda^{AB}_{22} \]

---

Constraint Sum-to-Zero

\[ \lambda^A_1 + \lambda^A_2 = 0 \]

\[ \lambda^B_1 + \lambda^B_2 = 0 \]

\[ \lambda^{AB}_{11} + \lambda^{AB}_{12} + \lambda^{AB}_{21} + \lambda^{AB}_{22} = 0 \]

---

Rumus Estimasi Parameter dengan Sum-to-Zero Constraint

\[ \lambda^A_1 = \frac{1}{2} \left[ (\log \mu_{11} + \log \mu_{12}) - (\log \mu_{21} + \log \mu_{22}) \right] \]

\[ \lambda^B_1 = \frac{1}{2} \left[ (\log \mu_{11} + \log \mu_{21}) - (\log \mu_{12} + \log \mu_{22}) \right] \]

\[ \lambda^{AB}_{12} = \frac{1}{4} \left[ \log \mu_{12} - \log \mu_{11} - \log \mu_{22} + \log \mu_{21} \right] \]

\[ \lambda_1^A = \frac{1}{2} \left[ (\log \mu_{11} + \log \mu_{12}) - (\log \mu_{21} + \log \mu_{22}) \right] \] \[ \lambda_1^B = \frac{1}{2} \left[ (\log \mu_{11} + \log \mu_{21}) - (\log \mu_{12} + \log \mu_{22}) \right] \] \[ \lambda_{12}^{AB} = \frac{1}{4} \left[ \log \mu_{12} - \log \mu_{11} - \log \mu_{22} + \log \mu_{21} \right] \] Analisis Data Tabel Kontingensi 2x2

Diberikan data:

Merokok	Sakit	Sehat
Ya	30	20
Tidak	10	40

Bentuk Model 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 \]

Definisi

• A1 = Merokok (Ya), A2 = Tidak
  
• B1 = Sakit, B2 = Sehat

---

Observasi

\[ \begin{aligned} n_{11} &= 30,\quad n_{12} = 20 \\ n_{21} &= 10,\quad n_{22} = 40 \end{aligned} \]

---

Sistem Persamaan Model Log-Linear

\[ \log(\mu_{11}) = \lambda + \lambda^A_1 + \lambda^B_1 + \lambda^{AB}_{11} \]

\[ \log(\mu_{12}) = \lambda + \lambda^A_1 + \lambda^B_2 + \lambda^{AB}_{12} \]

\[ \log(\mu_{21}) = \lambda + \lambda^A_2 + \lambda^B_1 + \lambda^{AB}_{21} \]

\[ \log(\mu_{22}) = \lambda + \lambda^A_2 + \lambda^B_2 + \lambda^{AB}_{22} \]

---

Constraint Sum-to-Zero

\[ \lambda^A_1 + \lambda^A_2 = 0 \]

\[ \lambda^B_1 + \lambda^B_2 = 0 \]

Langkah-langkah:

1. Hitung rata-rata log frekuensi sel:

\[ \lambda = \frac{1}{4} \sum_{i=1}^{2} \sum_{j=1}^{2} \log(n_{ij}) \]

\[ = \frac{1}{4} \left[ \log(30) + \log(20) + \log(10) + \log(40) \right] \]

\[ = 3.0971 \]

2. Efek utama A (Merokok):

\[ \lambda_1^A = \frac{1}{2} \left[ (\log(30) + \log(20)) - (\log(10) + \log(40)) \right] \]

\[ = \frac{1}{2} \left[ (3.4012 + 2.9957) - (2.3026 + 3.6889) \right] \]

\[ = \frac{1}{2} \left( 6.3969 - 5.9915 \right) \]

\[ = \frac{1}{2} \left( 0.4054 \right) \]

\[ = 0.2027 \]

\[ \lambda_2^A = -0.2027 \]

3. Efek utama B (Status):

\[ \lambda_{1}^{B} = \frac{1}{2}\left[(\log(30)+\log(10))-(\log(20)+\log(40))\right] \]

\[ = \frac{1}{2}\left[(3.4012+2.3026)-(2.9957+3.6889)\right] \]

\[ = \frac{1}{2}(5.7038-6.6846) \]

\[ = \frac{1}{2}(-0.9808) \]

\[ = -0.4904 \]

\[ \lambda_{2}^{B} = +0.4904 \]

4. Efek interaksi:

\[ \lambda_{11}^{AB} = \frac{1}{4}\left[\log(30)-\log(20)-\log(10)+\log(40)\right] \]

\[ = \frac{1}{4}\left[3.4012-2.9957-2.3026+3.6889\right] \]

\[ = \frac{1}{4}(3.4012-2.9957-2.3026+3.6889) \]

\[ = \frac{1}{4}(1.7918) \]

\[ = 0.4479 \]

\[ \lambda_{12}^{AB} = -\lambda_{11}^{AB}=-0.4479 \]

\[ \lambda_{21}^{AB} = -0.4479 \]

\[ \lambda_{22}^{AB} = +0.4479 \]

Ringkasan Parameter

• \(\lambda = 3.0971\)
• \(\lambda_1^A = 0.2027, \lambda_2^A = -0.2027\)
• \(\lambda_1^B = -0.4904, \lambda_2^B = 0.4904\)
• \(\lambda_{11}^{AB} = 0.4479, \lambda_{12}^{AB} = -0.4479, \lambda_{21}^{AB} = -0.4479, \lambda_{22}^{AB} = 0.4479\)

Hitung Odds Ratio dan Interval Kepercayaan

\[ OR = \frac{n_{11} n_{22}}{n_{12} n_{21}} = \frac{30 \times 40}{20 \times 10} = \frac{1200}{200} = 6 \] Log odds ratio: \[ \log(OR) = \log(6) = 1.7918 \] Standard error (SE):

\[ SE = \sqrt{ \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} } \\ \\ = \sqrt{ \frac{1}{30} + \frac{1}{20} + \frac{1}{10} + \frac{1}{40} }\\ \\ = \sqrt{0.0333 + 0.05 + 0.1 + 0.025}\\ \\ = \sqrt{0.2083}\\ \\ = 0.4564 \\ \]

95% Confidence Interval for \(\log(OR)\):

\[ \log(OR) \pm 1.96 \times SE = 1.7918 \pm 1.96 \times 0.4564 \]

\[ = (1.7918 - 0.895), \, (1.7918 + 0.895) \]

\[ = (0.8968, 2.6868) \]

Back-transform to get CI for \(OR\):

Lower = \(\exp(0.8968) = 2.452\)

Upper = \(\exp(2.6868) = 14.68\)

Jadi, OR=6 dengan taraf signifikansi 95% memiliki interval kepercayaan dalam rentang 2.452 - 14.68

Fitting Model Log-Linear dengan R

# Data 2x2
tabel <- matrix(c(30, 20, 10, 40), nrow = 2, byrow = TRUE)
colnames(tabel) <- c("Sakit", "Sehat")
rownames(tabel) <- c("Ya", "Tidak")
tabel

##       Sakit Sehat
## Ya       30    20
## Tidak    10    40

#Ubah menjadi bentuk tabel untuk glm
data <- as.data.frame(as.table(tabel))
colnames(data) <- c("Merokok", "Status", "Freq")
data

# Model tanpa interaksi
fit_no_inter <- glm(Freq ~ Merokok + Status, family = poisson, data = data)
summary(fit_no_inter)

## 
## Call:
## glm(formula = Freq ~ Merokok + Status, family = poisson, data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.996e+00  1.871e-01  16.013   <2e-16 ***
## MerokokTidak 3.892e-10  2.000e-01   0.000    1.000    
## StatusSehat  4.055e-01  2.041e-01   1.986    0.047 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 21.288  on 3  degrees of freedom
## Residual deviance: 17.261  on 1  degrees of freedom
## AIC: 43.036
## 
## Number of Fisher Scoring iterations: 4

# Model dengan interaksi
fit_inter <- glm(Freq ~ Merokok * Status, family = poisson, data = data)
summary(fit_inter)

## 
## Call:
## glm(formula = Freq ~ Merokok * Status, family = poisson, data = data)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.4012     0.1826  18.629  < 2e-16 ***
## MerokokTidak              -1.0986     0.3651  -3.009  0.00262 ** 
## StatusSehat               -0.4055     0.2887  -1.405  0.16015    
## MerokokTidak:StatusSehat   1.7918     0.4564   3.926 8.65e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.1288e+01  on 3  degrees of freedom
## Residual deviance: 3.9968e-15  on 0  degrees of freedom
## AIC: 27.775
## 
## Number of Fisher Scoring iterations: 3

Interpretasi Parameter

• Parameter Utama (Intercept) menunjukkan rata-rata log frekuensi sel
• Efek “Merokok” dan “Status”menunjukkan perbedaan log frekuensi antar kategori
• Interaksi signifika menunjukkan adanya asosiasi antara efek Merokok dan Status

Nilai \(\log(6)=1.7918\) sama dengan efek interaksi pada output R

Analisis Data Tabel Kontingensi 2x3

Suatu survei dilakukan untuk mengetahui hubungan antara Jenis Kelamin (Laki-laki/Perempuan) dan Kategori BMI (Kurus/Normal/Gemuk):

	Kurus	Normal	Gemuk
Laki-laki	12	20	8
Perempuan	18	24	10

Bentuk Model Log-Linear untuk Tabel 2x3

Bentuk umum model log-linear untuk tabel 2x3 (dengan sum-to-zero constraint): \[ \log(\mu_{ij}) = \lambda + \lambda_i^{A} + \lambda_j^{B} + \lambda_{ij}^{AB} \] dengan:

• \(\mu_{ij}\): ekspektasi frekuensi pada baris ke-i, kolom ke-j

- A: Jenis Kelamin (i=1:Laki-laki, i=2:Perempuan)

- B: Kategori BMI (j=1:Kurus, j=2:Normal, j=3:Gemuk)

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

Secara eksplisit :

\[ \log(\mu_{ij}) = \lambda \]

\[ + \lambda_1^A \, (\text{Laki-laki}), \, \lambda_2^A \, (\text{Perempuan}) \]

\[ + \lambda_1^B \, (\text{Kurus}), \, \lambda_2^B \, (\text{Normal}), \, \lambda_3^B \, (\text{Gemuk}) \]

\[ +\lambda_{ij}, (\text{interaksi jika ada}) \]

Fitting Model Log-Linear di R

# Membuat data frame dari tabel
tabel2x3 <- matrix(c(12, 20, 8, 18, 24, 10), nrow = 2, byrow = TRUE)
colnames(tabel2x3) <- c("Kurus", "Normal", "Gemuk")
rownames(tabel2x3) <- c("Laki-laki", "Perempuan")
tabel2x3

##           Kurus Normal Gemuk
## Laki-laki    12     20     8
## Perempuan    18     24    10

# Ubah menjadi data.frame untuk glm
data2x3 <- as.data.frame(as.table(tabel2x3))
colnames(data2x3) <- c("JenisKelamin", "BMI", "Freq")
data2x3

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

## 
## Call:
## glm(formula = Freq ~ JenisKelamin + BMI, family = poisson, data = data2x3)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             2.5683     0.2179  11.789   <2e-16 ***
## JenisKelaminPerempuan   0.2624     0.2103   1.248   0.2122    
## BMINormal               0.3830     0.2368   1.618   0.1058    
## BMIGemuk               -0.5108     0.2981  -1.713   0.0866 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 13.06443  on 5  degrees of freedom
## Residual deviance:  0.22527  on 2  degrees of freedom
## AIC: 35.26
## 
## Number of Fisher Scoring iterations: 3

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

## 
## Call:
## glm(formula = Freq ~ JenisKelamin * BMI, family = poisson, data = data2x3)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                       2.4849     0.2887   8.608   <2e-16 ***
## JenisKelaminPerempuan             0.4055     0.3727   1.088    0.277    
## BMINormal                         0.5108     0.3651   1.399    0.162    
## BMIGemuk                         -0.4055     0.4564  -0.888    0.374    
## JenisKelaminPerempuan:BMINormal  -0.2231     0.4802  -0.465    0.642    
## JenisKelaminPerempuan:BMIGemuk   -0.1823     0.6032  -0.302    0.762    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  1.3064e+01  on 5  degrees of freedom
## Residual deviance: -9.0719e-30  on 0  degrees of freedom
## AIC: 39.034
## 
## Number of Fisher Scoring iterations: 3

Interpretasi: Model tanpa interaksi: - Jika deviance tidak signifikan, maka Jenis Kelamin dan BMI independen.

• Intercept: log frekuensi pada kategori referensi (Laki-laki, Kurus)
• Koefisien JenisKelaminPerempuan: perbedaan log-frekuensi antara Perempuan vs Laki-laki (pada Kurus)
• Koefisien BMI: perbedaan log-frekuensi kategori BMI (Normal/Gemuk) terhadap Kurus (pada Laki-laki)

Model dengan interaksi:

• Jika koefisien interaksi signifikan, berarti ada hubungan/asosiasi antara Jenis Kelamin dan BMI. Artinya distribusi BMI berbeda antara Laki-laki dan Perempuan.

Contoh interpretasi hasil (misal): - Jika koefisien JenisKelaminPerempuan negatif: proporsi Perempuan pada kategori referensi lebih kecil dibanding Laki-laki. - Jika koefisien BMI_Normal positif: kemungkinan seseorang Normal lebih tinggi daripada Kurus (pada Laki-laki). - Jika model interaksi signifikan, pola distribusi BMI berbeda antara Laki-laki dan Perempuan.

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:

1. Model Saturated \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk} \] Model ini memuat semua kemungkinan interaksi, termasuk interaksi tiga arah (X, Y, dan Z).

2. Model Homogen \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \]

   Model ini hanya mengakomodasi interaksi dua arah antar variabel tanpa memasukkan interaksi tiga arah.

3. Model Conditional

   • Conditional pada X: \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} \] Memuat interaksi X dengan Y dan X dengan Z.

   • Conditional pada Y: \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk} \] Memuat interaksi Y dengan X dan Y dengan Z.

   • Conditional pada Z: \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} \] Memuat interaksi Z dengan X dan Z dengan Y.

4. Model Joint Independence

   • Independensi antara X & Y: \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XY}_{ij} \]

   • Independensi antara X & Z: \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XZ}_{ik} \]

   • Independensi antara Y & Z: \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{YZ}_{jk} \]

5. Model Tanpa Interaksi \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k \] Model ini hanya memasukkan efek utama tanpa interaksi antar variabel.

Pengujian Interaksi dalam Model Log-Linear 3 Arah

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

1. Pengujian Interaksi Tiga Arah (XYZ):
   • Bandingkan model saturated dengan model homogenous.
2. Pengujian Interaksi Dua Arah (XY, XZ, YZ):

   • Bandingkan model homogenous dengan model conditional.

   • Bandingkan model conditional dengan model joint independence.

   • Bandingkan model joint independence dengan model tanpa interaksi.

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

Soal Praktikum

Tabel berikut menyajikan data dari survei General Social Survey (GSS) tahun 1994 mengenai jenis kelamin responden, tingkat fundamentalisme, dan sikap terhadap hukuman mati untuk kasus pembunuhan. Susun dan interpretasikan model log-linear paling sederhana (paling parsimonious) untuk data ini. Jelaskan proses yang Anda lakukan dalam menentukan model terbaik serta asosiasi apa saja yang teridentifikasi. Tunjukkan juga bagaimana nilai yang diprediksi dari model menggambarkan asosiasi tersebut.

Fundamentalism	Jenis Kelamin	Mendukung	Menolak	Total
Fundamentalist	Laki-laki	128	32	160
Fundamentalist	Perempuan	123	73	196
Fundamentalist	Total	251	105	356
Moderate	Laki-laki	182	56	238
Moderate	Perempuan	168	105	273
Moderate	Total	350	161	511
Liberal	Laki-laki	119	49	168
Liberal	Perempuan	111	70	181
Liberal	Total	230	119	349

Keterangan: - Fundamentalisme: Fundamentalist, Moderate, Liberal - Jenis Kelamin: Laki-laki, Perempuan - Sikap: Mendukung (Favor), Menolak (Oppose) hukuman mati

Analisis Log-Linear untuk Tabel Tiga Arah

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

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

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

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

)
data

Membentuk Tabel Kontingensi 3 Arah

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

##                               Sikap 1fav 2opp
## Fundamentalisme Jenis_Kelamin                
## 1fund           1M                   128   32
##                 2F                   123   73
## 2mod            1M                   182   56
##                 2F                   168  105
## 3lib            1M                   119   49
##                 2F                   111   70

Analisis Log-Linear: Tahap Pemodelan Kita akan memodelkan tabel ini menggunakan beberapa model log-linear dan membandingkan kecocokan model (parsimonious model):

Uji Model Interaksi Tiga Arah (Saturated VS Homogenous)

Penentuan Kategori Referensi

##=============================##
# Penentuan kategori reference
##=============================##
x.sex  <- relevel(x.sex, ref = "2F")
y.fav  <- relevel(y.fav, ref = "2opp")
z.fund <- relevel(z.fund, ref = "3lib")

Model Saturated memasukkan semua interaksi hingga tiga arah: \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XY}_{ij} + \lambda^{XZ}_{ik} + \lambda^{YZ}_{jk} + \lambda^{XYZ}_{ijk} \]

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

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund + y.fav * z.fund + x.sex * y.fav * z.fund, 
##     family = poisson(link = "log"))
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    4.248495   0.119523  35.545  < 2e-16 ***
## x.sex1M                       -0.356675   0.186263  -1.915  0.05551 .  
## y.fav1fav                      0.461035   0.152626   3.021  0.00252 ** 
## z.fund1fund                    0.041964   0.167285   0.251  0.80193    
## z.fund2mod                     0.405465   0.154303   2.628  0.00860 ** 
## x.sex1M:y.fav1fav              0.426268   0.228268   1.867  0.06185 .  
## x.sex1M:z.fund1fund           -0.468049   0.282210  -1.659  0.09721 .  
## x.sex1M:z.fund2mod            -0.271934   0.249148  -1.091  0.27507    
## y.fav1fav:z.fund1fund          0.060690   0.212423   0.286  0.77511    
## y.fav1fav:z.fund2mod           0.008969   0.196903   0.046  0.96367    
## x.sex1M:y.fav1fav:z.fund1fund  0.438301   0.336151   1.304  0.19227    
## x.sex1M:y.fav1fav:z.fund2mod   0.282383   0.301553   0.936  0.34905    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.4536e+02  on 11  degrees of freedom
## Residual deviance: 5.9952e-15  on  0  degrees of freedom
## AIC: 100.14
## 
## Number of Fisher Scoring iterations: 3

exp(model_saturated$coefficients)

##                   (Intercept)                       x.sex1M 
##                    70.0000000                     0.7000000 
##                     y.fav1fav                   z.fund1fund 
##                     1.5857143                     1.0428571 
##                    z.fund2mod             x.sex1M:y.fav1fav 
##                     1.5000000                     1.5315315 
##           x.sex1M:z.fund1fund            x.sex1M:z.fund2mod 
##                     0.6262231                     0.7619048 
##         y.fav1fav:z.fund1fund          y.fav1fav:z.fund2mod 
##                     1.0625694                     1.0090090 
## x.sex1M:y.fav1fav:z.fund1fund  x.sex1M:y.fav1fav:z.fund2mod 
##                     1.5500717                     1.3262868

Model yang digunakan adalah model log-linear saturated dengan semua efek utama, interaksi dua arah, dan interaksi tiga arah. Model ini memodelkan hubungan antara jenis kelamin (x.sex), sikap terhadap hukuman mati (y.fav), dan tingkat fundamentalisme (z.fund) terhadap frekuensi responden. Dengan hasil estimasi koefisien nya sebagai berikut:

Parameter	Estimate	Std. Error	z value	Pr(>z)
(Intercept)	4.25	0.12	35.55	<2e-16***
x.sex1M	-0.36	0.19	-1.92	0.055 .
y.fav1fav	0.46	0.15	3.02	0.0025 **
z.fund1fund	0.04	0.17	0.25	0.80
z.fund2mod	0.41	0.15	2.63	0.0086 **
x.sex1M:y.fav1fav	0.43	0.23	1.87	0.062 .
x.sex1M:z.fund1fund	-0.47	0.28	-1.66	0.097 .
x.sex1M:z.fund2mod	-0.27	0.25	-1.09	0.28
y.fav1fav:z.fund1fund	0.06	0.21	0.29	0.78
y.fav1fav:z.fund2mod	0.01	0.20	0.05	0.96
x.sex1M:y.fav1fav:z.fund1fund	0.44	0.34	1.30	0.19
x.sex1M:y.fav1fav:z.fund2mod	0.28	0.30	0.94	0.35

Interpretasi Koefisien:

• (Intercept): Rata-rata log jumlah kasus untuk kategori referensi (Perempuan, Menolak hukuman mati, Liberal) adalah 4.25 (atau μ≈70 ).

• x.sex1M: Laki-laki memiliki expected count sekitar 0.7 kali Perempuan dalam kategori referensi lainnya, namun hanya mendekati signifikansi (p = 0.055).

• y.fav1fav: Mereka yang mendukung hukuman mati memiliki expected count sekitar 1.59 kali lipat dibanding yang menolak (signifikan, p = 0.0025).

• z.fund1fund: Kelompok Fundamentalist tidak berbeda nyata dari Liberal (exp(0.04)≈1.04; p = 0.80).

• z.fund2mod: Kelompok Moderate memiliki expected count 1.5 kali lebih besar dibanding Liberal (signifikan, p = 0.0086).

• Interaksi dua & tiga arah: Sebagian besar tidak signifikan (p > 0.05), artinya tidak ada bukti kuat adanya efek gabungan antar variabel.

Goodness-of-Fit:

• Residual deviance ≈0 menandakan model saturated benar-benar fit terhadap data (seluruh variasi data dijelaskan oleh model).

• AIC = 100.14 dapat digunakan untuk perbandingan dengan model yang lebih sederhana.

Kesimpulan:

• Model saturated ini sangat fit dengan data, namun tidak semua parameter/interaksi signifikan.

• Efek utama yang paling signifikan adalah:

  • Sikap mendukung hukuman mati (expected count 1.6x lebih tinggi dari yang menolak)

  • Kelompok Moderate (expected count 1.5x lebih tinggi dari Liberal)

  • Tidak ditemukan bukti kuat interaksi dua atau tiga arah yang signifikan.

• Model yang lebih sederhana (tanpa interaksi tiga arah) perlu dipertimbangkan untuk model final yang lebih parsimonious.

Catatan Interpretasi:

• Nilai exp(coef) menyatakan rasio ekspektasi (expected count ratio) dibandingkan baseline.

• Efek positif -> menaikkan expected count; Efek negatif -> menurunkan expected count.

• Koefisien signifikan pada p-value < 0.05

Model Homogenous

Model log-linear homogenous memasukkan semua efek utama dan semua interaksi dua arah, tanpa interaksi tiga arah. Secara matematis, model ini dapat dituliskan sebagai berikut: \[ \log(\mu_{ijk}) = \lambda + \lambda^{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.fav + z.fund +
              x.sex*y.fav + x.sex*z.fund + y.fav*z.fund,
              family = poisson(link = "log"))
summary(model_homogenous)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund + y.fav * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.31096    0.10522  40.972  < 2e-16 ***
## x.sex1M               -0.51575    0.13814  -3.733 0.000189 ***
## y.fav1fav              0.35707    0.12658   2.821 0.004788 ** 
## z.fund1fund           -0.06762    0.14452  -0.468 0.639854    
## z.fund2mod             0.33196    0.13142   2.526 0.011540 *  
## x.sex1M:y.fav1fav      0.66406    0.12728   5.217 1.81e-07 ***
## x.sex1M:z.fund1fund   -0.16201    0.15300  -1.059 0.289649    
## x.sex1M:z.fund2mod    -0.08146    0.14079  -0.579 0.562887    
## y.fav1fav:z.fund1fund  0.23873    0.16402   1.455 0.145551    
## y.fav1fav:z.fund2mod   0.13081    0.14951   0.875 0.381614    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.361  on 11  degrees of freedom
## Residual deviance:   1.798  on  2  degrees of freedom
## AIC: 97.934
## 
## Number of Fisher Scoring iterations: 3

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

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

Langkah-Langkah Pengujian: 1. Hipotesis: - H0: Tidak ada interaksi tiga arah (model homogenous sudah cukup) - H1: Ada interaksi tiga arah (model saturated diperlukan)

2. Hitung Selisih Deviance

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

## [1] 1.797977

3. Hitung Derajat Bebas

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

## [1] 2

4. Chi-Square Tabel (\(\alpha\)=0.05)

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

## [1] 5.991465

5. Keputusan Uji

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

## [1] "Terima H0"

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

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

Rangkuman

Pengujian ada tidaknya interaksi tiga arah (Saturated Model vs Homogenous Model)

• Hipotesis

  • H0: \(\lambda_{ijk}^{XYZ}=0\) (Tidak ada interaksi tiga arah; model yang terbentuk adalah model homogenous) -

  • H0: \(\lambda_{ijk}^{XYZ}\neq0\) (Tidak ada interaksi tiga arah; model yang terbentuk adalah model homogenous)

• Tingkat Signifikansi

  \(\alpha=5\%\)

• Statistik Uji

  • \(\Delta\)Deviance = Deviance model homogenous - Deviance model saturated = 1.798 - 0.00 = 1.798

  • \(db\) = \(db\) model homogenous - \(db\) model saturated = 2 - 0 = 2

• Daerah Penolakan

  Tolak H0 jika \(\Delta\)Deviance \(\>>\) \(\chi_{0.05,db}^2 = \chi_{0.05,2}^2 = 5.991\)

• Keputusan

  Karena 1.798 < 5.991, maka terima H0

• Interpretasi

  Pada taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi tiga arah antara jenis kelamin, fundamentalisme, dan pendapat mengenai hukuman mati.

Uji Model Interaksi Dua Arah (Homogenous vs 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.fav + z.fund +
              x.sex*y.fav + x.sex*z.fund,
              family = poisson(link = "log"))
summary(model_conditional_X)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          4.23495    0.08955  47.293  < 2e-16 ***
## x.sex1M             -0.52960    0.13966  -3.792 0.000149 ***
## y.fav1fav            0.48302    0.08075   5.982 2.20e-09 ***
## z.fund1fund          0.07962    0.10309   0.772 0.439916    
## z.fund2mod           0.41097    0.09585   4.288 1.81e-05 ***
## x.sex1M:y.fav1fav    0.65845    0.12708   5.181 2.20e-07 ***
## x.sex1M:z.fund1fund -0.12841    0.15109  -0.850 0.395405    
## x.sex1M:z.fund2mod  -0.06267    0.13908  -0.451 0.652274    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   3.9303  on  4  degrees of freedom
## AIC: 96.067
## 
## Number of Fisher Scoring iterations: 4

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

• Hipotesis

  • H0: \(\lambda_{jk}^{YZ}=0\) (Tidak ada interaksi antara pendapat hukuman mati (Y) dan fundamentalisme (Z)) -

  • H0: \(\lambda_{jk}^{YZ}\neq0\) (Ada interaksi antara pendapat hukuman mati (Y) dan fundamentalisme (Z))

• Tingkat Signifikansi

  \(\alpha=5\%\)

• Statistik Uji

  • \(\Delta\)Deviance = Deviance model conditional on X - Deviance model homogenous = 3.903 - 1.798 = 2.132

  • \(db\) = \(db\) model homogenous - \(db\) model saturated = 4 - 2 = 2

• Daerah Penolakan

  Tolak H0 jika \(\Delta\)Deviance \(\>>\) \(\chi_{0.05,db}^2 = \chi_{0.05,2}^2 = 5.991\)

• Keputusan

  Karena 2.132 < 5.991, maka terima H0

• Interpretasi

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

1. Hitung Selisih Deviance

# Pengujian hipotesis

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

## [1] 2.132302

2. Hitung Derajat Bebas

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

## [1] 2

3. Chi-Square Tabel

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

## [1] 5.991465

4. Keputusan Uji

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

## [1] "Terima"

Interpretasi Karena nilai Deviance.model = 2.13 lebih kecil dari nilai kritis chi-square tabel = 5.99 (dengan df = 2, alpha = 0.05), maka keputusan uji adalah “Terima”.

Pada taraf nyata 5%, belum cukup bukti untuk menolak H0, atau dengan kata lain tidak ada interaksi antara pendapat mengenai hukuman mati (Y) dan fundamentalisme (Z). Model yang terbentuk cukup hanya sampai dua interaksi dengan X (conditional on X), sehingga interaksi Y*Z tidak signifikan secara statistik.

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

Model log-linear conditional pada Y memasukkan efek utama dan interaksi dua arah antara X dengan Y dan Y dengan Z, tanpa interaksi antara X dengan Z maupun interaksi tiga arah. \[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XY}_{ij} + \lambda^{YZ}_{jk} \]

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

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     y.fav * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.33931    0.09919  43.748  < 2e-16 ***
## x.sex1M               -0.59345    0.10645  -5.575 2.48e-08 ***
## y.fav1fav              0.37259    0.12438   2.996  0.00274 ** 
## z.fund1fund           -0.12516    0.13389  -0.935  0.34989    
## z.fund2mod             0.30228    0.12089   2.500  0.01240 *  
## x.sex1M:y.fav1fav      0.65845    0.12708   5.181 2.20e-07 ***
## y.fav1fav:z.fund1fund  0.21254    0.16205   1.312  0.18966    
## y.fav1fav:z.fund2mod   0.11757    0.14771   0.796  0.42606    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   2.9203  on  4  degrees of freedom
## AIC: 95.057
## 
## Number of Fisher Scoring iterations: 4

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

• Hipotesis

  • H0: \(\lambda_{ik}^{XZ}=0\) (Tidak ada interaksi antara jenis kelamin (X) dan fundamentalisme (Z))

  • H0: \(\lambda_{ik}^{XZ}\neq0\) (Ada interaksi antara jenis kelamin

    24. dan fundamentalisme (Z))

• Tingkat Signifikansi

  \(\alpha=5\%\)

• Statistik Uji

  • \(\Delta\)Deviance = Deviance model conditional on Y - Deviance model homogenous = 2.9203 - 1.798 = 1.1223

  • \(db\) = \(db\) model homogenous - \(db\) model saturated = 4 - 2 = 2

• Daerah Penolakan

  Tolak H0 jika \(\Delta\)Deviance \(\>>\) \(\chi_{0.05,db}^2 = \chi_{0.05,2}^2 = 5.991\)

• Keputusan

  Karena 1.1223 < 5.991, maka terima H0

• Interpretasi

  Dengan taraf nyata 5%, belum cukup bukti untuk menolak H0 atau dapat dikatakan bahwa tidak ada interaksi antara jenis kelamin dan fundamentalisme. Dengan kata lain, model yang terbentuk adalah model tanpa parameter \(\lambda_{ik}^{XZ}\)

1. Hitung Selisih Deviance

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

## [1] 1.122315

2. Hitung Derajat Bebas

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

## [1] 2

3. Chi-Square Tabel

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

## [1] 5.991465

4. Keputusan Uji

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

## [1] "Terima"

Interpretasi Karena nilai Deviance.model = 1.12 lebih kecil dari nilai kritis chi-square tabel = 5.99 (dengan df = 2, alpha = 0.05), maka keputusan uji adalah “Terima”.

Pada taraf nyata 5%, belum cukup bukti untuk menolak H0, atau dengan kata lain tidak ada interaksi antara jenis kelamin (X) dan fundamentalisme (Z). Model yang terbentuk cukup 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 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.fav + z.fund +
              x.sex*z.fund + y.fav*z.fund,
              family = poisson(link = "log"))
summary(model_conditional_Z)

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * z.fund + 
##     y.fav * z.fund, family = poisson(link = "log"))
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            4.12255    0.10518  39.195  < 2e-16 ***
## x.sex1M               -0.07453    0.10713  -0.696    0.487    
## y.fav1fav              0.65896    0.11292   5.836 5.36e-09 ***
## z.fund1fund           -0.06540    0.15126  -0.432    0.665    
## z.fund2mod             0.33196    0.13777   2.410    0.016 *  
## x.sex1M:z.fund1fund   -0.12841    0.15109  -0.850    0.395    
## x.sex1M:z.fund2mod    -0.06267    0.13908  -0.451    0.652    
## y.fav1fav:z.fund1fund  0.21254    0.16205   1.312    0.190    
## y.fav1fav:z.fund2mod   0.11757    0.14771   0.796    0.426    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.361  on 11  degrees of freedom
## Residual deviance:  29.729  on  3  degrees of freedom
## AIC: 123.87
## 
## Number of Fisher Scoring iterations: 4

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

• Hipotesis

  • H0: \(\lambda_{ij}^{XY}=0\) (Tidak ada interaksi antara jenis kelamin (X) dan pendapat tentang hukuman mati (Y)) -

  • H0: \(\lambda_{ij}^{XY}\neq0\) (Ada interaksi antara jenis kelamin

    24. dan pendapat tentang hukuman mati (Y))

• Tingkat Signifikansi

  \(\alpha=5\%\)

• Statistik Uji

  • \(\Delta\)Deviance = Deviance model conditional on Y - Deviance model homogenous = 29.729 - 1.798 = 27.931

  • \(db\) = \(db\) model homogenous - \(db\) model saturated = 3 - 2 = 1

• Daerah Penolakan

  Tolak H0 jika \(\Delta\)Deviance \(\>>\) \(\chi_{0.05,db}^2 = \chi_{0.05,1}^2 = 3.841\)

• Keputusan

  Karena 27.931 > 3.841, maka tolak H0

• Interpretasi

  Dengan taraf nyata 5%, menolak H0 atau dapat dikatakan bahwa ada interaksi antara jenis kelamin dan pendapat tentang hukuman mati. Dengan kata lain, model yang terbentuk adalah model parameter \(\lambda_{ij}^{XY}\)

1. Hitung Selisih Deviance

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

## [1] 27.93095

2. Hitung Derajat Bebas

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

## [1] 1

3. Chi-Square Tabel

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

## [1] 3.841459

4. Keputusan Uji

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

## [1] "Tolak"

Interpretasi Karena nilai Deviance.model = 27.93 jauh lebih besar dari nilai kritis chi-square tabel = 3.84 (dengan df = 1, alpha = 0.05), maka keputusan uji adalah “Tolak”.

Pada taraf nyata 5%, terdapat cukup bukti untuk menolak H0, atau dengan kata lain ada interaksi antara jenis kelamin (X) dan pendapat tentang hukuman mati (Y). Model terbaik yang terbentuk adalah model dengan menyertakan parameter interaksi (conditional on Z).

Pemilihan Model Terbaik

Ringkasan Model Log Linear

Model	Parameter	Deviance	Jumlah Parameter	df	AIC
Saturated	λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢⱼˣʸ + λᵢₖˣᶻ + λⱼₖʸᶻ + λᵢⱼₖˣʸᶻ	0.000	12	0	100.14
Homogenous	λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢⱼˣʸ + λᵢₖˣᶻ + λⱼₖʸᶻ	1.798	10	2	97.934
Conditional on X	λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢⱼˣʸ + λᵢₖˣᶻ	3.9303	8	4	96.067
Conditional on Y	λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢⱼˣʸ + λⱼₖʸᶻ	2.9203	8	4	95.057
Conditional on Z	λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢₖˣᶻ + λⱼₖʸᶻ	29.729	9	3	123.87

Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah

Interaksi	Pengujian	\(\Delta\)Deviance	\(\Delta\)df	Chi-Square Tabel	Keputusan	Keterangan
XYZ	Saturated vs Homogenous	1.798	2	5.991	Terima H0	tidak ada interaksi
YZ	Conditional on X vs Homogenous	2.1323	2	5.991	Terima H0	tidak ada interaksi
XZ	Conditional on Y vs Homogenous	1.1223	2	5.991	Terima H0	tidak ada interaksi
XY	Conditional on Z vs Homogenous	27.931	1	3.841	Tolak H0	ada interaksi

Kesimpulan Pemilihan Model Terbaik

Dari hasil di atas diketahui bahwa asosiasi yang nyata hanya terdapat antara jenis kelamin dan pendapat mengenai hukuman mati. Sehingga, model terbaik adalah:

\[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XY}_{ij} \]

Model terbaik adalah model log-linear tanpa interaksi tiga arah dan hanya memuat interaksi dua arah antara jenis kelamin (X) dan sikap terhadap hukuman mati (Y).

Model Terbaik

Model terbaik dipilih berdasarkan pengujian interaksi yang signifikan, yaitu hanya interaksi dua arah antara jenis kelamin (X) dan sikap terhadap hukuman mati (Y):

\[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{XY}_{ij} \]

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

## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav, 
##     family = poisson(link = "log"))
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        4.26518    0.07794  54.721  < 2e-16 ***
## x.sex1M           -0.59345    0.10645  -5.575 2.48e-08 ***
## y.fav1fav          0.48302    0.08075   5.982 2.20e-09 ***
## z.fund1fund        0.01986    0.07533   0.264    0.792    
## z.fund2mod         0.38130    0.06944   5.491 4.00e-08 ***
## x.sex1M:y.fav1fav  0.65845    0.12708   5.181 2.20e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   4.6532  on  6  degrees of freedom
## AIC: 92.79
## 
## Number of Fisher Scoring iterations: 4

Model	AIC
Saturated	100.14
Homogenous	97.934
Conditional on X	96.067
Conditional on Y	95.057
Conditional on Z	123.87
Model Terbaik	92.79

Dari hasil di atas, terlihat bahwa model terbaik memiliki AIC yang lebih rendah dibandingkan dengan model lainnya seperti saturated, homogennous, dan conditional model.

Interpretasi Koefisien Model Terbaik

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

Koefisien	Nilai Koefisien	Exp(Nilai Koefisien)
Intercept	4.265	71.1776
x.sex1M	-0.5934	0.5524
y.fav1fav	0.483	1.621
z.fund1fund	0.0198	1.02
z.fund2mod	0.3813	1.464
x.sex1M:y.fav1fav	0.685	1.9318

• \(\exp(\lambda_{1M}^X) = \exp(-0.593) = 0.552 \rightarrow nilai odds\)

  Tanpa memperhatikan fundamentalisme dan pendapat mengenai hukuman mati, peluang seseorang berjenis kelamin laki-laki adalah 0,55 kali dibandingkan perempuan.
  Atau, peluang seseorang berjenis kelamin perempuan adalah \(1/0.55 = 1.81\) kali dibandingkan laki-laki.

• \(\exp(\lambda_{1fav}^Y) = \exp(0.483) = 1.621 \rightarrow nilai odds\)

  Tanpa memperhatikan jenis kelamin dan fundamentalisme, peluang seseorang mendukung hukuman mati adalah 1,621 kali dibandingkan yang menolak.

• \(\exp(\lambda_{1fund}^Z) = \exp(0.01986) = 1.02 \rightarrow nilai odds\)

  Tanpa memperhatikan jenis kelamin dan pendapat mengenai hukuman mati, peluang seseorang fundamentalist adalah 1,02 kali dibandingkan liberal.

• \(\exp(\lambda_{2mod}^Z) = \exp(0.381) = 1.464 \rightarrow nilai odds\)

  Tanpa memperhatikan jenis kelamin dan pendapat mengenai hukuman mati, peluang seseorang moderate adalah 1,464 kali dibandingkan liberal.

• \(\exp(\lambda_{1M,1fav}^XY) = \exp(0.658) = 1.932 \rightarrow nilai odds ratio\)

  Tanpa memperhatikan fundamentalisme, odds mendukung hukuman mati (dibandingkan menolak) jika dia laki-laki adalah 1,932 kali dibandingkan odds yang sama jika dia perempuan.

Nilai Dugaan Model Terbaik

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

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

Perhitungan Manual Nilai Dugaan (Fitted Value) Model Terbaik

\[ \hat{\mu}_{111} = \exp(\lambda + \lambda_{1m}^{x} + \lambda_{1fav}^{y} + \lambda_{fund}^{z} + \lambda_{1m,1fav}^{xy}) \] \[ = \exp(4.265 - 0.593 + 0.483 + 0.01986 + 0.658) \] \[ = \exp(4.833) = 125.595 \]

\[ \hat{\mu}_{112} = \exp(\lambda + \lambda_{1m}^{x} + \lambda_{1fav}^{y} + \lambda_{2mod}^{z} + \lambda_{1m,1fav}^{xy}) \] \[ = \exp(4.265 - 0.593 + 0.483 + 0.381 + 0.658) \] \[ = \exp(5.195) = 180.279 \]

\[ \hat{\mu}_{113} = \exp(\lambda + \lambda_{1m}^{x} + \lambda_{1fav}^{y} + \lambda_{lib}^{z} + \lambda_{1m,1fav}^{xy}) \] \[ = \exp(4.265 - 0.593 + 0.483 + 0 + 0.658) \] \[ = \exp(4.813) = 123.126 \]

\[ \hat{\mu}_{121} = \exp(\lambda + \lambda_{1m}^{x} + \lambda_{2opp}^{y} + \lambda_{fund}^{z} + \lambda_{1m,2opp}^{xy}) \] \[ = \exp(4.265 - 0.593 + 0 + 0.01986 + 0) \] \[ = \exp(3.692) = 40.109 \]

\[ \hat{\mu}_{122} = \exp(\lambda + \lambda_{1m}^{x} + \lambda_{2opp}^{y} + \lambda_{2mod}^{z} + \lambda_{1m,2opp}^{xy}) \] \[ = \exp(4.265 - 0.593 + 0 + 0.381 + 0) \] \[ = \exp(4.053) = 57.572 \]

\[ \hat{\mu}_{123} = \exp(\lambda + \lambda_{1m}^{x} + \lambda_{2opp}^{y} + \lambda_{lib}^{z} + \lambda_{1m,2opp}^{xy}) \] \[ = \exp(4.265 - 0.593 + 0 + 0 + 0) \] \[ = \exp(3.672) = 39.320 \]

\[ \hat{\mu}_{211} = \exp(\lambda + \lambda_{2f}^{x} + \lambda_{1fav}^{y} + \lambda_{fund}^{z} + \lambda_{2f,1fav}^{xy}) \] \[ = \exp(4.265 + 0 + 0.483 + 0.01986 + 0) \] \[ = \exp(4.768) = 117.691 \]

\[ \hat{\mu}_{212} = \exp(\lambda + \lambda_{2f}^{x} + \lambda_{1fav}^{y} + \lambda_{2mod}^{z} + \lambda_{2f,1fav}^{xy}) \] \[ = \exp(4.265 + 0 + 0.483 + 0.381 + 0) \] \[ = \exp(5.1295) = 168.933 \]

\[ \hat{\mu}_{213} = \exp(\lambda + \lambda_{2f}^{x} + \lambda_{1fav}^{y} + \lambda_{lib}^{z} + \lambda_{2f,1fav}^{xy}) \] \[ = \exp(4.265 + 0 + 0.483 + 0 + 0) \] \[ = \exp(4.748) = 115.377 \]

\[ \hat{\mu}_{221} = \exp(\lambda + \lambda_{2f}^{x} + \lambda_{2opp}^{y} + \lambda_{fund}^{z} + \lambda_{2f,2opp}^{xy}) \] \[ = \exp(4.265 + 0 + 0 + 0.01986 + 0) \] \[ = \exp(4.285) = 72.605 \]

\[ \hat{\mu}_{222} = \exp(\lambda + \lambda_{2f}^{x} + \lambda_{2opp}^{y} + \lambda_{2mod}^{z} + \lambda_{2f,2opp}^{xy}) \] \[ = \exp(4.265 + 0 + 0 + 0.381 + 0) \] \[ = \exp(4.646) = 104.217 \]

\[ \hat{\mu}_{223} = \exp(\lambda + \lambda_{2f}^{x} + \lambda_{2opp}^{y} + \lambda_{lib}^{z} + \lambda_{2f,2opp}^{xy}) \] \[ = \exp(4.265 + 0 + 0 + 0 + 0) \] \[ = \exp(4.265) = 71.178 \]

Keterangan: Nilai \(\hat{\mu}_{ijk}\) akan sama apapun referensi dari kategori peubahnya yang digunakan.

Referensi

Agresti, A. (2019). Categorical Data Analysis (3rd ed.). Wiley.

Fienberg, S. E. (2007). The Analysis of Cross-Classified Categorical Data. Springer.

Bishop, Y. M. M. (1975). Discrete Multivariate Analysis. MIT Press.

Gopinath, G., & Kumar, S. (2023). Understanding customer churn in streaming platforms: The role of user engagement and content consumption. Journal of Business Research, 157, 113576.

Hosmer, D. W., Jr., Lemeshow, S., & Sturdivant, R. X. (2013). Applied logistic regression (3rd ed.). John Wiley & Sons.

Dobson, A. J., & Barnett, A. G. (2018). An introduction to generalized linear models (4th ed.). Chapman & Hall/CRC.

Kim, M. S., & Kim, M. (2021). What drives consumption of video streaming services? A customer journey perspective. Electronic Commerce Research and Applications, 49, 101073.

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

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

Modul ADK, Pertemuan 13. (2025). Universitas Padjadjaran.