Tabel Kontingensi 2x2

Dalam analisis data, tabel kontingensi digunakan untuk memahami hubungan antara dua variabel kategorikal. Tabel kontingensi adalah tabel yang menyajikan distribusi frekuensi dari dua variabel dalam bentuk matriks. Tabel ini digunakan untuk menganalisis hubungan antara variabel dan menguji signifikansi statistik hubungan tersebut.

Tabel kontingensi 2x2 adalah bentuk paling sederhana dari tabel kontingensi, di mana terdapat dua kategori untuk masing-masing variabel. Dengan tabel ini, kita dapat menghitung berbagai ukuran probabilitas dan asosiasi antara variabel yang diamati.

Pada contoh ini, kita akan menganalisis hubungan antara kebiasaan olahraga dan kondisi kesehatan menggunakan tabel kontingensi 2x2.

Bentuk Umum Tabel Kontingensi 2x2:

Tabel kontingensi 2x2 biasanya disusun dalam format berikut:

\[ \begin{array}{|c|c|c|c|} \hline & \textbf{Kategori 1} & \textbf{Kategori 2} & \textbf{Total} \\ \hline \textbf{Kelompok 1} & n_{11} & n_{12} & n_{1.} \\ \textbf{Kelompok 2} & n_{21} & n_{22} & n_{2.} \\ \hline \textbf{Total} & n_{.1} & n_{.2} & n \\ \hline \end{array} \]

Dimana:

Distribusi Peluang Tabel Kontingensi 2x2

Dalam tabel kontingensi 2x2 terdapat beberapa distribusi peluang, yaitu peluang bersama, peluang marginal, dan peluang bersyarat.

Peluang Bersama

Peluang bersama adalah probabilitas bahwa dua kejadian terjadi secara bersamaan. Dalam tabel kontingensi 2x2, peluang bersama dihitung sebagai:

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

Dimana:

  • \(n_{ij}\) adalah jumlah observasi pada sel (i, j).

  • \(n\) adalah total jumlah observasi.

Peluang Marginal

Peluang marginal adalah probabilitas kejadian tunggal tanpa memperhatikan kejadian lain. Peluang marginal untuk suatu baris atau kolom dihitung sebagai:

\[ P(A_i) = \frac{n_{i+}}{n}, \quad P(B_j) = \frac{n_{+j}}{n} \]

Dimana:

  • \(n_{i+}\) adalah total observasi pada baris i.

  • \(n_{+j}\) adalah total observasi pada kolom j.

  • \(n\) adalah total jumlah observasi.

Peluang Bersyarat

Peluang bersyarat adalah probabilitas suatu kejadian terjadi dengan syarat kejadian lain sudah terjadi. Dalam tabel kontingensi, peluang bersyarat dihitung sebagai:

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

atau

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

Contoh Kasus

Sebuah tim peneliti dari sebuah universitas kesehatan masyarakat melakukan survei kepada 100 orang dewasa di kota X. Survei ini bertujuan untuk mengetahui apakah terdapat hubungan antara kebiasaan berolahraga dengan kondisi kesehatan seseorang. Dalam survei ini, responden diminta untuk mengisi kuesioner yang mencakup frekuensi mereka dalam melakukan aktivitas olahraga dan penilaian umum terhadap kondisi kesehatan mereka.

Peneliti mengelompokkan kebiasaan olahraga responden menjadi dua kategori: - Berolahraga - Tidak Berolahraga

Sementara itu, kondisi kesehatan dikelompokkan berdasarkan penilaian subyektif responden terhadap kesehatannya menjadi dua kategori: - Sehat - Tidak Sehat

Hasil dari survei ini disajikan dalam bentuk tabel kontingensi berikut: \[ \begin{array}{|c|c|c|c|} \hline & \textbf{Sehat (B1)} & \textbf{Tidak Sehat (B2)} & \textbf{Total} \\ \hline \textbf{Olahraga (A1)} & 30 & 10 & 40 \\ \textbf{Tidak Olahraga (A2)} & 20 & 40 & 60 \\ \hline \textbf{Total} & 50 & 50 & 100 \\ \hline \end{array} \]

  1. Hitung Peluang Bersama

    • \(P(A_1 \cap B_1) = \frac{30}{100} = 0.3\)

    • \(P(A_1 \cap B_2) = \frac{10}{100} = 0.1\)

    • \(P(A_2 \cap B_1) = \frac{20}{100} = 0.2\)

    • \(P(A_2 \cap B_2) = \frac{40}{100} = 0.4\)

Dengan menggunakan syntax R:

n_ij <- matrix(c(30, 10, 20, 40), nrow = 2, byrow = TRUE)
colnames(n_ij) <- c("Sehat", "Tidak Sehat")
rownames(n_ij) <- c("Olahraga", "Tidak Olahraga")
n <- sum(n_ij)
p_join <- n_ij / n
print(p_join)
##                Sehat Tidak Sehat
## Olahraga         0.3         0.1
## Tidak Olahraga   0.2         0.4
  1. Hitung Peluang Marginal

    • Peluang marginal untuk variabel olahraga:

      • \(P(A_1) = P(A_1 \cap B_1) + P(A_1 \cap B_2) = 0.3 + 0.1 = 0.4\)

      • \(P(A_2) = P(A_2 \cap B_1) + P(A_2 \cap B_2) = 0.2 + 0.4 = 0.6\)

    • Peluang marginal untuk variabel kesehatan:

      • \(P(B_1) = P(A_1 \cap B_1) + P(A_2 \cap B_1) = 0.3 + 0.2 = 0.5\)
      • \(P(B_2) = P(A_1 \cap B_2) + P(A_2 \cap B_2) = 0.1 + 0.4 = 0.5\)
n_ij <- matrix(c(30, 10, 20, 40), nrow = 2, byrow = TRUE)
p_marginal_A <- rowSums(p_join)
p_marginal_B <- colSums(p_join)
list(Peluang_Marginal_A = p_marginal_A, Peluang_Marginal_B = p_marginal_B)
## $Peluang_Marginal_A
##       Olahraga Tidak Olahraga 
##            0.4            0.6 
## 
## $Peluang_Marginal_B
##       Sehat Tidak Sehat 
##         0.5         0.5
  1. Hitung Peluang Bersyarat

    • \[ P(B1 | A1) = \frac{P(A1 \cap B1)}{P(A1)} = \frac{30}{40} = 0.75 \]

    • \[ P(B2 | A1) = \frac{P(A1 \cap B2)}{P(A1)} = \frac{10}{40} = 0.25 \]

    • \[ P(B1 | A2) = \frac{P(A2 \cap B1)}{P(A2)} = \frac{20}{60} = 0.33 \]

    • \[ P(B2 | A2) = \frac{P(A2 \cap B2)}{P(A2)} = \frac{40}{60} = 0.67 \]

    n_ij <- matrix(c(30, 10, 20, 40), nrow = 2, byrow = TRUE)
    colnames(n_ij) <- c("Sehat", "Tidak Sehat")
    rownames(n_ij) <- c("Olahraga", "Tidak Olahraga")
    p_conditional <- n_ij / rowSums(n_ij)
    print(p_conditional)
    ##                    Sehat Tidak Sehat
    ## Olahraga       0.7500000   0.2500000
    ## Tidak Olahraga 0.3333333   0.6666667

Interpretasi:

  • Peluang bersama menunjukkan peluang gabungan dari kejadian tertentu dalam tabel kontingensi 2x2

  • Peluang marginal menunjukkan peluang terjadinya suatu kejadian tanpa memperhatikan variabel lain

  • Peluang bersyarat menunjukkan peluang terjadinya suatu kejadian berubah ketika memperhatikan variabel lain

  • Dikarenakan, \(P(Sehat|Olahraga)=0.75 > P(Sehat|Tidak Olahraga)=0.25\), maka menunjukkan bahwa berolahraga lebih memberikan dampak sehat dibandingkan tidak berolahraga.

Ukuran Asosiasi Tabel Kontingensi 2x2

Ukuran asosiasi dalam tabel kontingensi 2x2 digunakan untuk mengukur hubungan antara dua variabel. Ukuran asosiasi yang umum digunakan, yaitu Risk Difference (RD)/Beda Peluang, Relative Risk/Resiko Relatif (RR), dan Odds Ratio (OR).

Risk Difference (RD)

Risk Difference (RD) atau beda peluang adalah selisih antara probabilitas kejadian dalam dua kelompok. RD menunjukkan seberapa besar perbedaan peluang antara kelompok yang terpapar dan kelompok yang tidak terpapar.

\[ RD = P(B1 | A1) - P(B1 | A2) \]

Jika RD = 0, berarti tidak terdapat perbedaan peluang antara dua kelompok tersebut.

Relative Risk (RR)

Relative Risk (RR) atau risiko relatif adalah perbandingan probabilitas kejadian dalam dua kelompok. RR digunakan untuk mengetahui seberapa besar risiko kejadian dalam kelompok terpapar dibandingkan dengan kelompok tidak terpapar.

\[ RR = \frac{P(B1 | A1)}{P(B1 | A2)} \]

Jika RR = 1, berarti tidak terdapat perbedaan risiko antara kedua kelompok tersebut.

Odds Ratio

Odds Ratio (OR) adalah ukuran yang membandingkan odds kejadian suatu peristiwa dalam dua kelompok. OR sering digunakan dalam studi kasus-kontrol.

\[ OR = \frac{(n_{11}/n_{12})}{(n_{21}/n_{22})} = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \] Dimana:

  • \(n_{11}\) = jumlah individu yang mengalami kejadian dalam kelompok pertama

  • \(n_{12}\) = jumlah individu yang tidak mengalami kejadian dalam kelompok pertama

  • \(n_{21}\) = jumlah individu yang mengalami kejadian dalam kelompok kedua

  • \(n_{22}\) = jumlah individu yang tidak mengalami kejadian dalam kelompok kedua

Jika OR = 1, tidak terdapat asosiasi antara 2 kelompok tersebut.

Contoh Kasus

Masih dengan kasus yang sama, kali ini akan dilihat lebih jauh seperti hubungan antara olahraga dengan kesehatan seseorang.

  1. Hitung Risk Difference (RD)

    \[ RD = \frac{30}{40} - \frac{20}{60} = 0.75 - 0.33 = 0.42 \]

    RD <- (n_ij[1,1] / sum(n_ij[1,])) - (n_ij[2,1] / sum(n_ij[2,]))
    print(RD)
    ## [1] 0.4166667
  2. Hitung Relative Risk (RR)

    \[ RR = \frac{30/40}{20/60} = \frac{0.75}{0.33} = 2.27 \]

    RR <- (n_ij[1,1] / sum(n_ij[1,])) / (n_ij[2,1] / sum(n_ij[2,]))
    print(RR)
    ## [1] 2.25
  3. Hitung Odds Ratio (OR)

    \[ OR = \frac{(30 \times 40)}{(10 \times 20)} = \frac{1200}{200} = 6 \]

    OR <- (n_ij[1,1] * n_ij[2,2]) / (n_ij[1,2] * n_ij[2,1])
    print(OR)
    ## [1] 6

Interpretasi:

  • RD = 0.42, maka orang yang berolahraga memiliki peluang 0.42 lebih tinggi untuk sehat dibandingkan dengan orang yang tidak berolahraga.

  • RR = 2.27, maka orang yang berolahraga 2.27 kali lebih mungkin untuk sehat dibandingkan dengan orang yang tidak berolahraga.

  • OR = 6, maka orang yang berolahraga 6 kali lebih mungkin untuk sehat dibandingkan dengan orang yang tidak berolahraga.

Inferensi Tabel Kontingensi Dua Arah

Dalam statistik, inferensi adalah suatu proses pengampilan kesimpulan mengenai populasi berdasarkan sampel data. Dalam konteks tabel kontingensi dua arah, inferensi dilakukan untuk menganalisis hubungan antara dua variabel kategorik yang ada dalam tabel kontingensi. Terdapat 2 jenis inferensi dalam tabel kontingensi, yaitu estimasi dan pengujian.

Estimasi

Estimasi statistik adalah proses untuk memperkirakan nilai suatu parameter populasi berdasarkan data dari sampel. Karena data populasi sering tidak tersedia secara keseluruhan, kita menggunakan data sampel untuk memperkirakan parameter populasi seperti rata-rata (mean), proporsi (proportion), atau variansi (variance).

Estimasi terbagi menjadi dua jenis utama:

  • Estimasi Titik (Point Estimation)
  • Estimasi Interval (Interval Estimation)

Estimasi Titik

Estimasi titik memberikan perkiraan tunggal terhadap suatu parameter populasi.

Contoh:

  • Estimasi rata-rata populasi \(\mu\) menggunakan rata-rata sampel \(\bar{x}\)
  • Estimasi proporsi populasi \(p\) menggunakan proporsi sampel \(\hat{p}\)

Notasi:

\[ \bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i \]

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

di mana: - \(\bar{x}\): rata-rata sampel - \(\hat{p}\): proporsi sampel - \(x_i\): data ke-i - \(n\): ukuran sampel - \(x\): jumlah kejadian “sukses”

Estimasi Interval

Estimasi interval memberikan rentang nilai di sekitar estimasi titik, yang diduga kuat mencakup nilai parameter populasi dengan tingkat kepercayaan tertentu, seperti 95%.

Interval Kepercayaan untuk Rata-rata (σ tidak diketahui)

Jika data berdistribusi normal dan simpangan baku populasi tidak diketahui, maka digunakan distribusi t-Student:

\[ \bar{x} \pm t_{\alpha/2, \, df} \cdot \frac{s}{\sqrt{n}} \]

Keterangan: - \(\bar{x}\): rata-rata sampel
- \(s\): simpangan baku sampel
- \(n\): ukuran sampel
- \(df = n - 1\): derajat bebas
- \(t_{\alpha/2}\): nilai t pada tingkat kepercayaan tertentu

Interval Kepercayaan untuk Proporsi

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

Keterangan: - \(\hat{p}\): proporsi sampel - \(z_{\alpha/2}\): nilai z dari distribusi normal standar - \(n\): ukuran sampel

Uji Hipotesis

Contoh Kasus

Sebuah kelompok peneliti di bidang kesehatan masyarakat ingin mengetahui apakah konsumsi vitamin C secara rutin dapat mengurangi risiko terkena flu. Untuk itu, mereka melakukan studi terhadap 200 orang dewasa selama musim flu berlangsung (sekitar 3 bulan).

Responden dibagi menjadi dua kelompok berdasarkan kebiasaan mereka mengonsumsi suplemen vitamin C:

  • Mengonsumsi Vitamin C secara rutin
  • Tidak mengonsumsi Vitamin C sama sekali

Setelah musim flu berakhir, para peneliti mencatat apakah masing-masing responden terkena flu atau tidak terkena flu selama periode pengamatan. Data hasil pengamatan disusun dalam tabel kontingensi 2x2 berikut:

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Terkena Flu (B1)} & \textbf{Tidak Terkena Flu (B2)} & \textbf{Total} \\ \hline \textbf{Minum Vitamin (A1)} & 25 & 75 & 100 \\ \textbf{Tidak Minum Vitamin (A2)} & 50 & 50 & 100 \\ \hline \textbf{Total} & 75 & 125 & 200 \\ \hline  \end{array} \]

Uji Proporsi

Uji proporsi dalam tabel kontingensi 2x2 digunakan untuk membandingkan proporsi dua kelompok dan menentukan apakah perbedaan antara proporsi tersebut signifikan secara statistik.

  • Hipotesis:

    • Hipotesis nol (H₀): Tidak ada perbedaan proporsi antara dua kelompok (\(p_1 = p_2\)).
    • Hipotesis alternatif (H₁): Ada perbedaan proporsi antara dua kelompok (\(p_1 \neq p_2\)).
  • Statistik Uji:

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

    • \(\hat{p}_1\) dan \(\hat{p}_2\) adalah proporsi sampel dari masing-masing kelompok.

    • \(p\) adalah proporsi gabungan dari kedua kelompok.

    • \(n_1\) dan \(n_2\) adalah ukuran sampel masing-masing kelompok.

    Dengan,

    • \(\hat{p}_1 = \frac{n_{11}}{n_{1.}}\)

    • \(\hat{p}_2 = \frac{n_{21}}{n_{2.}}\)

    • \(\hat{p} = \frac{n_{11}+n_{21}}{n_{1.}+n_{2.}}\)

  • Kriteria Uji:

    • Jika nilai Z lebih besar dari nilai kritis dalam distribusi normal standar, maka H₀ ditolak, yang berarti terdapat perbedaan signifikan antara dua proporsi.

    • Jika nilai p-value lebih kecil dari tingkat signifikansi (misalnya 0.05), maka H₀ ditolak.

  1. Hipotesis

    • Hipotesis nol (H₀): Tidak ada perbedaan proporsi antara dua kelompok (\(p_1 = p_2\)).
    • Hipotesis alternatif (H₁): Ada perbedaan proporsi antara dua kelompok (\(p_1 \neq p_2\)).
  2. Hitung Proporsi Sampel

    \[ \hat{p}_1 = \frac{25}{100} = 0.25, \quad \hat{p}_2 = \frac{50}{100} = 0.50 \]

  3. Hitung Proporsi Gabungan

    \[ \hat{p} = \frac{25 + 50}{100 + 100} = \frac{75}{200} = 0.375 \]

  4. Hitung Standar Error

    \[ SE = \sqrt{p(1-p) \left(\frac{1}{n_1} + \frac{1}{n_2}\right)} \] \[ SE = \sqrt{0.375 \times 0.625 \left(\frac{1}{100} + \frac{1}{100}\right)} \] \[ SE = \sqrt{0.2344 \times 0.02} = \sqrt{0.004688} = 0.0685 \]

  5. Hitung Statistik Uji

    \[ Z = \frac{\hat{p}_1 - \hat{p}_2}{SE} \] \[ Z = \frac{0.25 - 0.50}{0.0685} = \frac{-0.25}{0.0685} = -3.65 \]

  6. Tentukan Kriteria Uji

    Untuk tingkat signifikansi \(\alpha = 0.05\), nilai kritis dari distribusi normal standar adalah 1.96.

    Karena nilai Z = -3.65 lebih kecil dari -1.96, maka H₀ ditolak, yang berarti terdapat perbedaan signifikan antara kelompok yang mengonsumsi vitamin dan yang tidak dalam kemungkinan terkena flu.

Dengan menggunakan syntax R:

data_vitamin <- matrix(c(25, 75, 50, 50), nrow = 2, byrow = TRUE)
rownames(data_vitamin) <- c("Minum Vitamin", "Tidak Minum Vitamin")
colnames(data_vitamin) <- c("Terkena Flu", "Tidak Terkena Flu")
# Uji proporsi untuk data konsumsi vitamin
test_result_vitamin <- prop.test(x = c(25, 50), n = c(100, 100), alternative = "two.sided", correct = FALSE)

# Menampilkan hasil uji proporsi
test_result_vitamin
## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  c(25, 50) out of c(100, 100)
## X-squared = 13.333, df = 1, p-value = 0.0002607
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.3796394 -0.1203606
## sample estimates:
## prop 1 prop 2 
##   0.25   0.50

Interpretasi:

  • Dari hasil perhitungan di atas, didapatkan bahwa terdapat perbedaan signifikan antara kelompok yang mengonsumsi vitamin C dan yang tidak mengonsumsi vitamin C terhadap kemungkinan terkena flu.

  • Maka, konsumsi vitamin C memiliki pengaruh signifikan terhadap kemungkinan terkena flu. Seseorang yang mengonsumsi vtamin C memiliki kemungkinan lebih kecil terkena flu dibandingkan dengan yang tidak mengonsumsi vitamin C

Uji Asosiasi

Untuk menguji signifikansi dari ketiga ukuran asosiasi, kita dapat menggunakan statistik Z. Asumsinya adalah data berasal dari sampel besar sehingga pendekatan normal dapat digunakan. Dengan hipotesis sebagai berikut:

  • Hipotesis nol (H₀): Tidak ada asosiasi antara dua kelompok

  • Hipotesis alternatif (H₁): Ada asosiasi antara dua kelompok

Statistik Uji:

  1. Risk Difference (RD)

    • Standard error Risk Difference (RD)

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

    • Statistik uji Risk Difference (RD)

    \[ Z_{RD} = \frac{RD}{SE_{RD}} \]

  2. Relative Risk (RR)

    • Standard error Relative Risk (RR)

      \[ SE_{\log RR} = \sqrt{ \frac{1}{n_{11}} - \frac{1}{n_{1.}} + \frac{1}{n_{21}} - \frac{1}{n_{2.}} } \]

    • Statistik Uji Relative Risk (RR)

    \[ Z_{RR} = \frac{\log(RR)}{SE_{\log RR}} \]

  3. Odds Ratio (OR)

    • Standard error Odds Ratio (OR)

      \[ SE_{\log OR} = \sqrt{ \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} } \]

    • Statistik Uji Z Odds Ratio

      \[ Z_{OR} = \frac{\log(OR)}{SE_{\log OR}} \]

Kriteria uji:

  • Jika nilai Z lebih besar dari nilai kritis dalam distribusi normal standar, maka H₀ ditolak, yang berarti terdapat asosiasi anatara 2 kelompok.

  • Jika nilai p-value lebih kecil dari tingkat signifikansi (misalnya 0.05), maka H₀ ditolak.

Contoh kasus:

Misalkan dilakukan penelitian untuk mengetahui apakah konsumsi vitamin C dapat mengurangi kemungkinan seseorang terkena flu. Berikut adalah hasil penelitian dalam tabel kontingensi 2x2:

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Terkena Flu (B1)} & \textbf{Tidak Terkena Flu (B2)} & \textbf{Total} \\ \hline \textbf{Minum Vitamin (A1)} & 25 & 75 & 100 \\ \textbf{Tidak Minum Vitamin (A2)} & 50 & 50 & 100 \\ \hline \textbf{Total} & 75 & 125 & 200 \\ \hline  \end{array} \]

  1. Risk Difference (RD)

    \[ RD = \hat{p}_1 - \hat{p}_2 = 0.25 - 0.50 = -0.25 \]

    \[ SE_{RD} = \sqrt{ \frac{\hat{p}_1(1 - \hat{p}_1)}{n_{1.}} + \frac{\hat{p}_2(1 - \hat{p}_2)}{n_{2.}} } = \sqrt{ \frac{0.25 \times 0.75}{100} + \frac{0.50 \times 0.50}{100} } = \sqrt{ 0.001875 + 0.0025 } = \sqrt{0.004375} \approx 0.0661 \]

    \[ Z_{RD} = \frac{RD}{SE_{RD}} = \frac{-0.25}{0.0661} \approx -3.78 \]

  2. Relative Risk (RR)

    \[ RR = \frac{\hat{p}_1}{\hat{p}_2} = \frac{0.25}{0.50} = 0.5 \]

    \[ SE_{\log RR} = \sqrt{ \frac{1}{n_{11}} - \frac{1}{n_{1.}} + \frac{1}{n_{21}} - \frac{1}{n_{2.}} } = \sqrt{ \frac{1}{25} - \frac{1}{100} + \frac{1}{50} - \frac{1}{100} } = \sqrt{0.04 - 0.01 + 0.02 - 0.01} = \sqrt{0.04} = 0.2 \]

    \[ Z_{RR} = \frac{\log(RR)}{SE_{\log RR}} = \frac{\log(0.5)}{0.2} = \frac{-0.6931}{0.2} \approx -3.47 \]

  3. Odds Ratio (OR)

    \[ OR = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} = \frac{25 \cdot 50}{75 \cdot 50} = \frac{1250}{3750} = 0.333 \]

    \[ SE_{\log OR} = \sqrt{ \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} } = \sqrt{ \frac{1}{25} + \frac{1}{75} + \frac{1}{50} + \frac{1}{50} } = \sqrt{0.04 + 0.0133 + 0.02 + 0.02} = \sqrt{0.0933} \approx 0.3055 \]

    \[ Z_{OR} = \frac{\log(OR)}{SE_{\log OR}} = \frac{\log(0.333)}{0.3055} = \frac{-1.0986}{0.3055} \approx -3.60 \]

Dengan menggunakan syntax R:

# Data
a <- 25  # Minum vitamin, terkena flu
b <- 75  # Minum vitamin, tidak flu
c <- 50  # Tidak minum, terkena flu
d <- 50  # Tidak minum, tidak flu

n1 <- a + b  # Total yang minum
n2 <- c + d  # Total yang tidak minum

# Hitung proporsi risiko
p1 <- a / n1  # Risiko kelompok minum
p2 <- c / n2  # Risiko kelompok tidak minum

## 1. Risk Difference (RD)
RD <- p1 - p2
SE_RD <- sqrt( (p1 * (1 - p1)) / n1 + (p2 * (1 - p2)) / n2 )
Z_RD <- RD / SE_RD
p_RD <- 2 * pnorm(-abs(Z_RD))

## 2. Risk Ratio (RR)
RR <- p1 / p2
SE_logRR <- sqrt((1 - p1) / (a) + (1 - p2) / (c))
Z_logRR <- log(RR) / SE_logRR
p_RR <- 2 * pnorm(-abs(Z_logRR))

## 3. Odds Ratio (OR)
OR <- (a * d) / (b * c)
SE_logOR <- sqrt(1/a + 1/b + 1/c + 1/d)
Z_logOR <- log(OR) / SE_logOR
p_OR <- 2 * pnorm(-abs(Z_logOR))

## Tampilkan hasil
cat("=== RISK DIFFERENCE (RD) ===\n")
## === RISK DIFFERENCE (RD) ===
cat("RD:", RD, "\nZ:", Z_RD, "\np-value:", p_RD, "\n\n")
## RD: -0.25 
## Z: -3.779645 
## p-value: 0.0001570523
cat("=== RISK RATIO (RR) ===\n")
## === RISK RATIO (RR) ===
cat("RR:", RR, "\nZ:", Z_logRR, "\np-value:", p_RR, "\n\n")
## RR: 0.5 
## Z: -3.465736 
## p-value: 0.0005287824
cat("=== ODDS RATIO (OR) ===\n")
## === ODDS RATIO (OR) ===
cat("OR:", OR, "\nZ:", Z_logOR, "\np-value:", p_OR, "\n")
## OR: 0.3333333 
## Z: -3.596053 
## p-value: 0.0003230822

Interpretasi:

  • Untuk tingkat signifikansi \(\alpha = 0.05\), nilai kritis dari distribusi normal standar adalah 1.96.

    Karena nilai Z = -3.78, -3.47, -3.60 lebih kecil dari -1.96, maka H₀ ditolak, yang berarti terdapat asosiasi antara 2 kelompok.

  • Maka, antara kelompok konsumsi vitamin C memiliki asosiasi kelompok kemungkinan terkena flu.

Uji Independensi

Uji independensi pada tabel kontingensi 2x2 adalah metode statistik untuk menentukan apakah dua variabel kategori saling bebas (independen) atau berhubungan (dependen) satu sama lain. Tujuannya untuk menentukan apakah ada asosiasi atau ketergantungan antara variabel baris (A) dan kolom (B).

Uji Chi-Square

  1. Hipotesis

    • Hipotesis Nol (H₀): Kedua variabel independen (tidak saling berhubungan).

    • Hipotesis Alternatif (H₁): Kedua variabel tidak independen (ada hubungan/asosiasi).

  2. Statistik Uji:

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

    Dengan:

    • \(O\): frekuensi pengamatan (observed)
    • \(E\): frekuensi harapan (expected)

    Frekuensi harapan dihitung sebagai:

    \[ E_{ij} = \frac{(\text{total baris}_i) \cdot (\text{total kolom}_j)}{n} \]

  3. Kriteria Uji:

    • Jika \(\chi^2_{\text{hitung}} > \chi^2_{\text{tabel}}\), maka tolak H₀

    • Jika \(\chi^2_{\text{hitung}} \leq \chi^2_{\text{tabel}}\), maka gagal tolak H₀

    • Jika p-value < 0.05 → Tolak \(H_0\), artinya ada hubungan antara dua variabel.

    • Jika p-value ≥ 0.05 → Gagal tolak \(H_0\), artinya tidak ada bukti cukup bahwa kedua variabel saling berhubungan.

Contoh Kasus:

Misalkan dilakukan penelitian untuk mengetahui apakah konsumsi vitamin C dapat mengurangi kemungkinan seseorang terkena flu. Berikut adalah hasil penelitian dalam tabel kontingensi 2x2:

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Terkena Flu (B1)} & \textbf{Tidak Terkena Flu (B2)} & \textbf{Total} \\ \hline \textbf{Minum Vitamin (A1)} & 25 & 75 & 100 \\ \textbf{Tidak Minum Vitamin (A2)} & 50 & 50 & 100 \\ \hline \textbf{Total} & 75 & 125 & 200 \\ \hline  \end{array} \]

  1. Hipotesis:

    • H₀ (Hipotesis Nol): Tidak ada hubungan antara minum vitamin dan kejadian flu (variabel independen).

    • H₁ (Hipotesis Alternatif): Ada hubungan antara minum vitamin dan kejadian flu (variabel tidak independen).

  2. Statistik uji:

    • Hitung Nilai Frekuensi Harapan:

      \[ E_{ij} = \frac{(\text{Total Baris}_i) \times (\text{Total Kolom}_j)}{n} \]

      • \(E_{11} = \frac{100 \times 75}{200} = 37.5\)
      • \(E_{12} = \frac{100 \times 125}{200} = 62.5\)
      • \(E_{21} = \frac{100 \times 75}{200} = 37.5\)
      • \(E_{22} = \frac{100 \times 125}{200} = 62.5\)
    • Hitung Statistik uji Chi-Square

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

      \[\begin{align*} \chi^2 &= \frac{(25 - 37.5)^2}{37.5} + \frac{(75 - 62.5)^2}{62.5} + \frac{(50 - 37.5)^2}{37.5} + \frac{(50 - 62.5)^2}{62.5} \\ &= \frac{(-12.5)^2}{37.5} + \frac{(12.5)^2}{62.5} + \frac{(12.5)^2}{37.5} + \frac{(-12.5)^2}{62.5} \\ &= \frac{156.25}{37.5} + \frac{156.25}{62.5} + \frac{156.25}{37.5} + \frac{156.25}{62.5} \\ &= 4.1667 + 2.5 + 4.1667 + 2.5 \\ &= 13.3334 \end{align*}\]

    • Hitung nilai Chi-Square Tabel

      Derajat kebebasan: \(df = (2 - 1)(2 - 1) = 1\) - Taraf signifikansi: \(\alpha = 0.05\) - Nilai kritis dari tabel Chi-Square:
      \[ \chi^2_{(0.05;1)} = 3.841 \]

Dengan menggunakan syntax R:

# Buat tabel kontingensi
flu_matrix <- matrix(c(25, 75, 50, 50), nrow = 2, byrow = TRUE)

# Tambahkan nama baris dan kolom
colnames(flu_matrix) <- c("Flu", "Tidak_Flu")
rownames(flu_matrix) <- c("Minum_Vitamin", "Tidak_Minum")
flu_table <- as.table(flu_matrix)

# Uji chi-square independensi
chi_result <- chisq.test(flu_table, correct = FALSE)  # continuity correction dinonaktifkan untuk hasil lebih murni

# Tampilkan hasil
print(chi_result)
## 
##  Pearson's Chi-squared test
## 
## data:  flu_table
## X-squared = 13.333, df = 1, p-value = 0.0002607
  1. Keputusan

    • Untuk tingkat signifikansi \(\alpha = 0.05\), nilai kritis dari distribusi chi-square adalah 3.841.

      Karena nilai X2 = 13.3334 lebih besar dari 3.841, maka H₀ ditolak, yang berarti terdapat hubungan antara variabel konsumsi vitamin C dengan kemungkinan terkena flu.

  2. Kesimpulan

    Maka, antara kelompok konsumsi vitamin C memiliki hubungan kelompok kemungkinan terkena flu dan tidak saling bebas (tidak independen).

Uji Partisi Chi-Square

Partisi chi-square adalah prosedur penguraian (dekomposisi) nilai chi-square total dari tabel kontingensi berukuran I × J (baris × kolom) menjadi komponen-komponen chi-square yang lebih kecil. Ini berguna untuk mengetahui sumber utama penyimpangan antara data yang diamati dan yang diharapkan.

Misalnya: jika hasil uji chi-square menunjukkan adanya hubungan antara dua variabel, partisi chi-square dapat membantu menjelaskan dimana tepatnya hubungan tersebut terjadi, yaitu di baris atau kolom tertentu.

Langkah-langkah:

  • Hitung nilai Chi-Square total

  • Jika signifikan, maka lanjutkan ke partisi

  • Tentukan strategi partisi

  • Hitung nilai chi-square dengan data yang sudah di partisi

  • Bandingkan dengan nilai kritis.

Contoh Kasus:

Sebuah tim peneliti sosial dari sebuah universitas ingin mengetahui apakah terdapat hubungan antara jenis kelamin dan preferensi minuman favorit masyarakat. Penelitian ini dilakukan dengan menyebarkan kuesioner kepada 150 responden dari berbagai kalangan usia dan profesi.

Dalam kuesioner, responden diminta untuk memilih salah satu jenis minuman yang paling mereka sukai dari tiga pilihan: - Kopi - Teh - Jus

Selain itu, identitas jenis kelamin juga dicatat, dan dikelompokkan ke dalam dua kategori: - Pria - Wanita

Setelah data dikumpulkan dan diklasifikasikan, hasilnya ditampilkan dalam tabel kontingensi berikut:

\[\begin{array}{|c|c|c|c|c|} \hline \textbf{Jenis Kelamin} & \textbf{Kopi} & \textbf{Teh} & \textbf{Jus} & \textbf{Total} \\ \hline \text{Pria} & 30 & 20 & 10 & 60 \\ \text{Wanita} & 10 & 30 & 20 & 60 \\ \hline \text{Total} & 40 & 50 & 30 & 120 \\ \hline \end{array}\]

# Matriks Observasi
observed <- matrix(c(30, 20, 10, 10, 30, 20), 
                   nrow = 2, byrow = TRUE,
                   dimnames = list(
                     Jenis_Kelamin = c("Pria", "Wanita"),
                     Minuman = c("Kopi", "Teh", "Jus")
                   ))
observed
##              Minuman
## Jenis_Kelamin Kopi Teh Jus
##        Pria     30  20  10
##        Wanita   10  30  20
  1. Hitung nilai chi-square total

    \[ \begin{array}{|c|c|c|c|} \hline E_{ij} & \text{Kopi} & \text{Teh} & \text{Jus} \\ \hline \text{Pria} & \frac{60 \times 40}{120} = 20 & \frac{60 \times 50}{120} = 25 & \frac{60 \times 30}{120} = 15 \\ \text{Wanita} & 20 & 25 & 15 \\ \hline \end{array} \]

    \[\begin{align*} \chi^2 &= \frac{(30 - 20)^2}{20} + \frac{(20 - 25)^2}{25} + \frac{(10 - 15)^2}{15} \\ &\quad + \frac{(10 - 20)^2}{20} + \frac{(30 - 25)^2}{25} + \frac{(20 - 15)^2}{15} \\ &= \frac{100}{20} + \frac{25}{25} + \frac{25}{15} + \frac{100}{20} + \frac{25}{25} + \frac{25}{15} \\ &= 5 + 1 + 1.67 + 5 + 1 + 1.67 \\ &= 15.34 \end{align*}\]

    Dengan derajat bebas:

    \[ df = (2 - 1)(3 - 1) = 2 \]

    Bandingkan nilai \(\chi^2 = 15.34\) dengan nilai kritis \(\chi^2_{0.05, 2} = 5.991\)

    Karena \(15.34 > 5.991\), maka terdapat hubungan signifikan antara jenis kelamin dan preferensi minuman.

# Hitung frekuensi harapan
expected <- chisq.test(observed)$expected
expected
##              Minuman
## Jenis_Kelamin Kopi Teh Jus
##        Pria     20  25  15
##        Wanita   20  25  15
# Uji chi-square dari R
chi_test <- chisq.test(observed, correct = FALSE)
chi_test
## 
##  Pearson's Chi-squared test
## 
## data:  observed
## X-squared = 15.333, df = 2, p-value = 0.0004682
  • Melanjutkan ke partisi chi-square

    Buat tabel 2x2 baru:

    Kopi+Teh Jus Total
    Pria 50 10 60
    Wanita 40 20 60
    Total 90 30 120
    # Gabungkan kolom Kopi dan Teh menjadi satu
    observed_partisi <- matrix(c(50, 10, 40, 20), 
                           nrow = 2, byrow = TRUE,
                           dimnames = list(
                             Jenis_Kelamin = c("Pria", "Wanita"),
                             Minuman = c("Kopi+Teh", "Jus")
                           ))
    observed_partisi
    ##              Minuman
    ## Jenis_Kelamin Kopi+Teh Jus
    ##        Pria         50  10
    ##        Wanita       40  20

    Hitung nilai chi-square partisi:

    \[ E_{11} = \frac{60 \times 90}{120} = 45,\quad E_{12} = \frac{60 \times 30}{120} = 15 \]

    \[ \chi^2_{\text{partisi}} = \frac{(50 - 45)^2}{45} + \frac{(10 - 15)^2}{15} + \frac{(40 - 45)^2}{45} + \frac{(20 - 15)^2}{15} \]

    \[ = \frac{25}{45} + \frac{25}{15} + \frac{25}{45} + \frac{25}{15} = 0.56 + 1.67 + 0.56 + 1.67 = 4.46 \]

expected_partisi <- chisq.test(observed_partisi)$expected
expected_partisi
##              Minuman
## Jenis_Kelamin Kopi+Teh Jus
##        Pria         45  15
##        Wanita       45  15
chi_partisi_test <- chisq.test(observed_partisi, correct = FALSE)
chi_partisi_test
## 
##  Pearson's Chi-squared test
## 
## data:  observed_partisi
## X-squared = 4.4444, df = 1, p-value = 0.03501
  • Dengan \(df = 1\), nilai kritis \(\chi^2_{0.05,1} = 3.841\)

    Karena \(4.46 > 3.841\), maka partisi ini signifikan, artinya penyimpangan utama terjadi antara kelompok yang memilih Jus dan lainnya.

  • Interpretasi:

    • Uji chi-square menunjukkan hubungan signifikan.

    • Partisi chi-square mengindikasikan bahwa penyimpangan utama terjadi pada pilihan minuman Jus antara pria dan wanita.

Uji Likelihood Ratio (G2)

Uji Likelihood Ratio (G2) adalah alternatif dari uji Chi-Square Pearson yang digunakan untuk menguji apakah dua variabel kategori bersifat independen. Dalam konteks tabel kontingensi 2x2, uji ini berguna untuk memverifikasi apakah ada asosiasi antara dua variabel.

Hipotesis

  • \(H_0\): Baris dan kolom bersifat independen
  • \(H_1\): Baris dan kolom tidak independen

Rumus Statistik Likelihood Ratio

Rumus untuk menghitung nilai statistik Likelihood Ratio (G²) pada tabel kontingensi adalah sebagai berikut:

\[ G^2 = 2 \sum_{i=1}^{r} \sum_{j=1}^{c} O_{ij} \cdot \ln \left( \frac{O_{ij}}{E_{ij}} \right) \]

Di mana:

  • \(O_{ij}\): Frekuensi observasi pada sel ke-(i,j)
  • \(E_{ij}\): Frekuensi harapan pada sel ke-(i,j) yang dihitung dengan:

\[ E_{ij} = \frac{(\text{Total Baris}_i) \cdot (\text{Total Kolom}_j)}{N} \]

  • \(G^2\): Statistik likelihood ratio (juga disebut deviance)
  • \(\ln\): Logaritma natural
  • \(r\): Jumlah baris
  • \(c\): Jumlah kolom

Derajat Bebas

Untuk tabel kontingensi berukuran 2x2:

\[ df = (r - 1)(c - 1) = (2 - 1)(2 - 1) = 1 \]

Kriteria Uji:

  • Jika \(G^2_{\text{hitung}} > \chi^2_{\text{tabel}}\), maka tolak H₀

  • Jika \(G^2_{\text{hitung}} \leq \chi^2_{\text{tabel}}\), maka gagal tolak H₀

  • Jika p-value < 0.05 → Tolak \(H_0\), artinya ada hubungan antara dua variabel.

  • Jika p-value ≥ 0.05 → Gagal tolak \(H_0\), artinya tidak ada bukti cukup bahwa kedua variabel saling berhubungan.

Contoh Kasus:

Sebuah penelitian ingin mengetahui apakah terdapat hubungan antara jenis kelamin dan status merokok.

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Merokok (B1)} & \textbf{Tidak Merokok (B2)} & \textbf{Total} \\ \hline \textbf{Pria (A1)} & 30 & 20 & 50 \\ \textbf{Wanita (A2)} & 10 & 40 & 50 \\ \hline \textbf{Total} & 40 & 60 & 100 \\ \hline \end{array} \]

# Data Observasi
observed <- matrix(c(30, 20, 10, 40),
                   nrow = 2, byrow = TRUE,
                   dimnames = list(
                     Jenis_Kelamin = c("Pria", "Wanita"),
                     Merokok = c("Merokok", "Tidak_Merokok")
                   ))
observed
##              Merokok
## Jenis_Kelamin Merokok Tidak_Merokok
##        Pria        30            20
##        Wanita      10            40
# Marginal dan total
row_totals <- rowSums(observed)
col_totals <- colSums(observed)
grand_total <- sum(observed)

row_totals
##   Pria Wanita 
##     50     50
col_totals
##       Merokok Tidak_Merokok 
##            40            60
grand_total
## [1] 100
  1. Hitung frekuensi ekspektasi

    • \(E_{11} = \frac{50 \cdot 40}{100} = 20\)

    • \(E_{12} = \frac{50 \cdot 60}{100} = 30\)

    • \(E_{21} = \frac{50 \cdot 40}{100} = 20\)

    • \(E_{22} = \frac{50 \cdot 60}{100} = 30\)

# Hitung nilai E_ij manual
E11 <- (row_totals[1] * col_totals[1]) / grand_total
E12 <- (row_totals[1] * col_totals[2]) / grand_total
E21 <- (row_totals[2] * col_totals[1]) / grand_total
E22 <- (row_totals[2] * col_totals[2]) / grand_total

expected_manual <- matrix(c(E11, E12, E21, E22),
                          nrow = 2, byrow = TRUE,
                          dimnames = list(
                            Jenis_Kelamin = c("Pria", "Wanita"),
                            Merokok = c("Merokok", "Tidak_Merokok")
                          ))
expected_manual
##              Merokok
## Jenis_Kelamin Merokok Tidak_Merokok
##        Pria        20            30
##        Wanita      20            30
  1. Hitung statistik uji (G2)

    Hitung satu per satu:

    • \(30 \cdot \ln\left(\frac{30}{20}\right) = 30 \cdot \ln(1.5) = 30 \cdot 0.4055 = 12.165\)
    • \(20 \cdot \ln\left(\frac{20}{30}\right) = 20 \cdot \ln(0.6667) = 20 \cdot (-0.4055) = -8.11\)
    • \(10 \cdot \ln\left(\frac{10}{20}\right) = 10 \cdot \ln(0.5) = 10 \cdot (-0.6931) = -6.931\)
    • \(40 \cdot \ln\left(\frac{40}{30}\right) = 40 \cdot \ln(1.3333) = 40 \cdot 0.2877 = 11.508\)

    Jumlah total:

    \[ G^2 = 2 \cdot (12.165 - 8.11 - 6.931 + 11.508) = 2 \cdot 8.632 = 17.264 \]

# Hitung G² secara manual
G2_manual <- 2 * (
  observed[1,1] * log(observed[1,1] / expected_manual[1,1]) +
  observed[1,2] * log(observed[1,2] / expected_manual[1,2]) +
  observed[2,1] * log(observed[2,1] / expected_manual[2,1]) +
  observed[2,2] * log(observed[2,2] / expected_manual[2,2])
)
G2_manual
## [1] 17.26092
  1. Hitung derajat bebas

    \[ df = (r - 1)(c - 1) = (2 - 1)(2 - 1) = 1 \]

    p_value_manual <- pchisq(G2_manual, df = 1, lower.tail = FALSE)
    p_value_manual
    ## [1] 3.258188e-05
  2. Keputusan

    Nilai \(G^2 = 17.264\), dibandingkan dengan nilai kritis \(\chi^2_{0.05,1} = 3.841\). Karena:

    \[ G^2 = 17.264 > 3.841 \]

    Dan \[ \text{p-value} \approx 0.00003 \]

    Karena p-value < 0.05 atau \(G^2\) > nilai kritis, maka kita menolak \(H_0\).

  3. Kesimpulan

    Maka, antara kelompok jenis kelamin memiliki hubungan dengan kelompok merokok (tidak independen).

Uji Eksak Fisher

Uji Independensi Eksak Fisher (Fisher’s Exact Test) digunakan untuk menguji apakah dua variabel kategorik pada tabel kontingensi berukuran kecil (umumnya 2x2) saling bebas (independen) atau tidak. Uji ini sangat tepat digunakan ketika frekuensi harapan kecil, khususnya jika ada sel dengan nilai harapan kurang dari 5.

Uji ini tidak mengandalkan distribusi chi-square dan memberikan nilai p eksak, berbeda dari pendekatan aproksimasi.

Hipotesis:

  • H₀ (Hipotesis nol): Kedua variabel bersifat independen (tidak ada hubungan).

  • H₁ (Hipotesis alternatif): Kedua variabel tidak independen (ada hubungan).

Rumus Distribusi Hipergeometrik:

Untuk menghitung peluang konfigurasi tabel kontingensi tertentu, digunakan:

\[ P(X = n_{11}) = \frac{{\binom{n_{1.}}{n_{11}} \cdot \binom{n_{2.}}{n_{.1} - n_{11}}}}{{\binom{n}{n_{.1}}}} \]

Keterangan:

  • \(n_{1.} =\) jumlah total baris 1
  • \(n_{2.} =\) jumlah total baris 2
  • \(n_{.1} =\) jumlah total kolom 1
  • \(n =\) total observasi
  • \(n_{11} =\) sel [1,1]

Contoh Kasus:

Sebuah studi ingin melihat hubungan antara jenis kelamin dan preferensi minuman:

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Kopi (B1)} & \textbf{Teh (B2)} & \textbf{Total} \\ \hline \textbf{Pria (A1)} & 8 & 2 & 10 \\ \textbf{Wanita (A2)} & 1 & 5 & 6 \\ \hline \textbf{Total} & 40 & 60 & 100 \\ \hline \end{array} \]

Kita ingin menghitung peluang terjadinya tabel yang sama ekstrem atau lebih ekstrem daripada tabel observasi:

  • Peluang untuk \(X = 8\) (tabel observasi):

\[ P(X = 8) = \frac{\binom{10}{8} \cdot \binom{6}{1}}{\binom{16}{9}} = \frac{45 \cdot 6}{11440} = \frac{270}{11440} \approx 0.0236 \]

  • Peluang untuk \(X = 9\) (lebih ekstrem mendukung H₁):

    \[ P(X = 9) = \frac{\binom{10}{9} \cdot \binom{6}{0}}{\binom{16}{9}} = \frac{10 \cdot 1}{11440} = \frac{10}{11440} \approx 0.00087 \]

  • Peluang untuk \(X = 3\) (lebih ekstrem mendukung H₁):

    \[ P(X = 3) = \frac{\binom{10}{3} \cdot \binom{6}{6}}{\binom{16}{9}} = \frac{120 \cdot 1}{11440} = \frac{10}{11440} \approx 0.01049 \]

  • Total p-value (dua sisi konservatif):

    \[ p\text{-value} = P(X=8) + P(X=9) + P(X=3) = 0.0236 + 0.00087 + 0.01049 = \boxed{0.03496} \]

Jika, menggunakan fisher.test() dengan R:

# Data tabel kontingensi
data <- matrix(c(8, 2, 1, 5), 
               nrow = 2, 
               byrow = TRUE,
               dimnames = list("Jenis Kelamin" = c("Pria", "Wanita"),
                               "Minuman" = c("Kopi", "Teh")))

data
##              Minuman
## Jenis Kelamin Kopi Teh
##        Pria      8   2
##        Wanita    1   5
# Uji Fisher
fisher.test(data, alternative = "two.sided")
## 
##  Fisher's Exact Test for Count Data
## 
## data:  data
## p-value = 0.03497
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##     1.008849 1049.791446
## sample estimates:
## odds ratio 
##   15.46969

Kesimpulan:

Dengan p-value sebesar 0.03496 < 0.05, kita menolak H₀. Terdapat bukti bahwa terdapat hubungan antara jenis kelamin dan preferensi minuman dalam sampel tersebut.

Tabel Kontingensi 3 Arah

Tabel kontingensi 3 arah adalah ekstensi dari tabel kontingensi dua arah yang digunakan untuk menyajikan dan menganalisis hubungan antara tiga variabel kategorik. Dalam bentuk ini, data disusun dalam array tiga dimensi dengan ukuran \(I \times J \times K\), di mana:

Tujuan Analisis:

Analisis tabel kontingensi 3 arah bertujuan untuk:

Struktur Tabel Kontingensi 3 Arah:

Dengan tabel marjinal sebagai berikut:

Total frekuensi pada masing-masing variabel bisa dihitung dengan:

Contoh Kasus:

Data menampilkan tabel kontingensi 3 arah, yang memperlihatkan hubungan pengaruh jenis kelamin, konsumsi minuman beralkohol terhadap perilaku seks bebas pada remaja. Data berjumlah 620 dengan proporsi 350 laki-laki, dan 270 perempuan.

\[\begin{array}{|c|c|c|c|} \hline \textbf{Jenis Kelamin} & \textbf{Alkohol} & \textbf{Seks Bebas: Ya} & \textbf{Seks Bebas: Tidak} \\ \hline Laki-Laki & Ya & 130 & 100 \\ Laki-Laki & Tidak & 30 & 90 \\ Perempuan & Ya & 70 & 90 \\ Perempuan & Tidak & 20 & 90 \\ \hline \end{array}\]

Tabel Parsial dan Marginal

Tabel marginal: jumlah total untuk kombinasi dua variabel (X dan Y) dengan mengabaikan yang ketiga (Z). Tabel parsial: tabel dua arah pada kondisi tetap satu variabel (Z) tertentu (kondisional).

Tabel Frekuensi Parsial

Tabel frekuensi parsial adalah tabel dua arah yang menunjukkan hubungan antara dua variabel (X dan Y) dengan mengendalikan (memfiksasi) nilai dari variabel ketiga (Z).

Tabel Frekuensi parsial untuk Z (Jenis Kelamin) Laki-Laki:* \[\begin{array}{|c|c|c|c|} \hline & \textbf{Seks Bebas (Ya)} & \textbf{Seks Bebas (Tidak)} & \textbf{Total} \\ \hline \textbf{Alkohol} & 130 & 100 & 230 \\ \textbf{Tidak Alkohol} & 30 & 90 & 120 \\ \hline \textbf{Total} & 160 & 190 & 350 \\ \hline \end{array}\]

Tabel Frekuensi parsial untuk Z (Jenis Kelamin) Perempuan:

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Seks Bebas (Ya)} & \textbf{Seks Bebas (Tidak)} & \textbf{Total} \\ \hline \textbf{Alkohol} & 70 & 90 & 160 \\ \textbf{Tidak Alkohol} & 20 & 90 & 110 \\ \hline \textbf{Total} & 90 & 180 & 270 \\ \hline \end{array}\]

Dengan menggunakan syntax R:

data3 <- array(c(130, 30, 100, 90, 70, 20, 90, 90 ),
dim = c(2, 2, 2),
dimnames = list(
Alkohol = c("Ya", "Tidak"),
Seks_Bebas = c("Ya", "Tidak"),
Jenis_Kelamin = c("Laki-Laki", "Perempuan")
))
# Ekstrak tabel parsial berdasarkan jenis kelamin
freq_parsial_laki <- data3[, , "Laki-Laki"]
freq_parsial_perempuan <- data3[, , "Perempuan"]
# Tampilkan hasil
freq_parsial_laki
##        Seks_Bebas
## Alkohol  Ya Tidak
##   Ya    130   100
##   Tidak  30    90
freq_parsial_perempuan
##        Seks_Bebas
## Alkohol Ya Tidak
##   Ya    70    90
##   Tidak 20    90

Tabel frekuensi parsial adalah alat penting dalam analisis data kategorikal 3 arah, karena memungkinkan kita memahami dinamika hubungan antar variabel secara lebih terperinci dan kontekstual. Ini membantu menghindari generalisasi yang menyesatkan dari tabel total (agregat) dan mengarahkan pada interpretasi yang lebih akurat.

Tabel Frekuensi Marginal

Tabel frekuensi marginal adalah tabel yang menunjukkan jumlah total dari suatu variabel dengan mengabaikan (mengakumulasi) variabel lain dalam tabel kontingensi. Pada tabel 3 arah yang melibatkan tiga variabel kategorikal (misalnya X, Y, dan Z), frekuensi marginal diperoleh dengan menjumlahkan frekuensi atas satu atau dua variabel.

Tabel Frekuensi Marginal untuk X (Alkohol) dan Y (Seks Bebas): \[\begin{array}{|c|c|c|c|} \hline & \textbf{Seks Bebas} & \textbf{Tidak Seks Bebas} & \textbf{Total} \\ \hline \textbf{Alkohol} & 200 & 190 & 390 \\ \textbf{Tidak Alkohol} & 50 & 180 & 230 \\ \hline \textbf{Total} & 250 & 370 & 620 \\ \hline \end{array}\]

Tabel Frekuensi Marginal untuk Z (Jenis Kelamin) dan Y (Seks Bebas):

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Seks Bebas} & \textbf{Tidak Seks Bebas} & \textbf{Total} \\ \hline \textbf{Laki-Laki} & 160 & 190 & 350 \\ \textbf{Perempuan} & 90 & 180 & 270 \\ \hline \textbf{Total} & 350 & 270 & 620 \\ \hline \end{array}\]

Dengan menggunakan syntax R:

data4 <- array(c(130, 30, 100, 90, 70, 20, 90, 90 ),
dim = c(2, 2, 2),
dimnames = list(
Alkohol = c("Ya", "Tidak"),
Seks_Bebas = c("Ya", "Tidak"),
Jenis_Kelamin = c("Laki-Laki", "Perempuan")
))
# Ekstrak tabel marginal
freq_marginal_X <- apply(data4,1,sum)
freq_marginal_Z <- apply(data4,3,sum)

freq_marginal_X
##    Ya Tidak 
##   390   230
freq_marginal_Z
## Laki-Laki Perempuan 
##       350       270

Kesimpulan:

Tabel frekuensi marginal sangat penting dalam analisis data kategorikal karena memberikan informasi global tentang distribusi dua variabel. Memberikan gambaran umum distribusi dua variabel tanpa mempertimbangkan variabel ketiga. Membantu mengevaluasi pengaruh keseluruhan dari satu variabel terhadap yang lain.

Distribusi Peluang

Peluang Bersama

Peluang bersama menunjukkan kemungkinan terjadinya ketiga peristiwa secara bersamaan:

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

Sebagai contoh digunakan data pada bagian sebelumnya mengenai hubungan jenis kelamin, konsumsi alkohol, dan seks bebas. Peluang laki-laki, konumsi alkohol, dan perilaku seks bebas:

\[ P(Laki-Laki, Alkohol,Seks Bebas) = \frac{130}{620} = 0.20967 \]

Dengan menggunakan syntax R:

data <- data.frame(
  Jenis_Kelamin = c(rep("Laki-laki", 2), rep("Perempuan", 2)),
  Konsumsi_Alkohol = rep(c("Ya", "Tidak"), 2),
  Seks_Bebas_Ya = c(130, 30, 70, 20),
  Seks_Bebas_Tidak = c(100, 90, 90, 90),
  Total = c(230, 120, 160, 110)
)
data$P_ZXY <- data$Seks_Bebas_Ya/sum(data$Total)
data

Tabel Probabilitas Bersama diperoleh dengan membagi semua nilai di dalam sel dengan jumlah keseluruhan sampel.

Tabel Probabilitas Bersama untuk Y (Seks Bebas) dan Tidak Seks Bebas:

\[\begin{array}{|c|c|cc|} \hline \textbf{Z (Jenis Kelamin)}& \textbf{X (Alkohol)} & \textbf{P(Y} & \textbf{X,Z)}\\ \hline \textbf{Laki-Laki} & Ya & 0.20967 & 0.16129 \\ \textbf{Laki-Laki} & Tidak & 0.04839 & 0.14516 \\ \textbf{Perempuan} & Ya & 0.11290 & 0.03226 \\ \textbf{Perempuan} & Tidak & 0.14516 & 0.14516 \\ \hline \end{array}\]

Dengan menggunakan syntax R:

data5 <- array(c(130, 30, 100, 90, 70, 20, 90, 90 ),
dim = c(2, 2, 2),
dimnames = list(
Alkohol = c("Ya", "Tidak"),
Seks_Bebas = c("Ya", "Tidak"),
Jenis_Kelamin = c("Laki-Laki", "Perempuan")
))

##Hitung probabilitas bersama
total <- sum(data5)
p_bersama <- data5/total
ftable(p_bersama)
##                    Jenis_Kelamin  Laki-Laki  Perempuan
## Alkohol Seks_Bebas                                    
## Ya      Ya                       0.20967742 0.11290323
##         Tidak                    0.16129032 0.14516129
## Tidak   Ya                       0.04838710 0.03225806
##         Tidak                    0.14516129 0.14516129

Peluang Marginal

Peluang marginal menunjukkan peluang terjadinya satu atau dua peristiwa, tanpa mempertimbangkan variabel lainnya. Untuk menghitung peluang marginal dapat dengan menjumlahkan probabilitas bersama.

\[ P(\text{Seks Bebas} = \text{Ya}) = \frac{130+30+70+20}{620} = \frac{250}{620} \approx 0.403225 \]

\[ P(\text{Seks Bebas} = \text{Tidak}) = \frac{100+90+90+90}{620} = \frac{370}{620} \approx 0.59677 \]

Dengan menggunakan syntax R:

P_Y1 <- sum(data$Seks_Bebas_Ya)/sum(data$Total)
P_Y1
## [1] 0.4032258
P_Y2 <- sum(data$Seks_Bebas_Tidak)/sum(data$Total)
P_Y2
## [1] 0.5967742

Tabel Peluang Marginal dihitung dengan menjumlahkan probabilitas bersama untuk setiap variabel yang tersisa setelah mengabaikan variabel lain.

Tabel marginal X (Alkohol) dan Y (Seks Bebas) \[\begin{array}{|c|c|cc|} \hline \textbf{X (Alkohol)} & \textbf{P(Y)} & \textbf{P(Tidak Y)}\\ \hline \textbf{Ya} & 0.30645 & 0.32258 \\ \textbf{Tidak} & 0.08065 & 0.29032\\ \hline \end{array}\]

\[\begin{array}{|c|c|cc|} \hline \textbf{Z (Jenis Kelamin)}& \textbf{P(Y)} & \textbf{P(Tidak Y)}\\ \hline \textbf{Laki-Laki} & 0.25806 & 0.30645 \\ \textbf{Perempuan} & 0.14516 & 0.29032 \\ \hline \end{array}\]

Dengan menggunakan syntax R:

data6 <- array(c(130, 30, 100, 90, 70, 20, 90, 90 ),
dim = c(2, 2, 2),
dimnames = list(
Alkohol = c("Ya", "Tidak"),
Seks_Bebas = c("Ya", "Tidak"),
Jenis_Kelamin = c("Laki-Laki", "Perempuan")
))

# Hitung probabilitas Bersama
total <- sum(data6)
p_bersama <- data6/total

# Hitung probabilitas marginal
marginal_X <- apply(p_bersama,1,sum)
marginal_Z <- apply(p_bersama,3,sum)

marginal_X
##        Ya     Tidak 
## 0.6290323 0.3709677
marginal_Z
## Laki-Laki Perempuan 
## 0.5645161 0.4354839

Kesimpulan

Tabel peluang marginal membantu memahami distribusi probabilitas dari masing-masing variabel tanpa mempertimbangkan variabel lainnya. Ini memberikan wawasan penting dalam analisis data kategori dengan tabel kontingensi tiga arah.

Peluang Bersyarat

Peluang bersyarat didefinisikan sebagai: \[ P(Z|X,Y) = \frac{P(Z,X,Y)}{P(X,Y)} \]

\[\begin{align*} &\text{Total populasi: } N = 620 \\ \\ &\text{Probabilitas bersama: } \\ &P(Z = Laki-Laki, X = \text{Ya}, Y = \text{Ya}) = \frac{130}{620} \\ &P(X = \text{Ya}, Y = \text{Ya}) = \frac{130+70}{620} = \frac{200}{620} \\ \\ &\text{Peluang bersyarat:} \\ &P(Z = Laki-Laki \mid X = \text{Ya}, Y = \text{Ya}) = \frac{P(Z = Laki-Laki, X = \text{Ya}, Y = \text{Ya})}{P(X = \text{Ya}, Y = \text{Ya})} \\ &= \frac{130 / 620}{200 / 620} = \frac{130}{200} \approx 0.65 \\ \end{align*}\]

Dengan menggunakan syntax R:

P_X_Y1 <- sum(data$Seks_Bebas_Ya [data$Konsumsi_Alkohol == "Ya"]) / sum(data$Total)
P_Laki_given_X_Y1 <- (130 / 620) / P_X_Y1
P_Laki_given_X_Y1
## [1] 0.65

Kesimpulan:

  • Peluang bersama dihitung dengan membagi frekuensi dengan total populasi data

  • Peluang marginal dihitung dengan menjumlahkan probabilitas bersama

  • Peluang bersyarat dihitung dengan menggunakan rumus Bayes

Ukuran Asosiasi

Ukuran Asosiasi digunakan untuk mengukur kekuatan hubungan antara dua variabel kategori. Terdapat 3 ukuran asosiasi yang umum digunakan, yaitu:

  • Risk Difference

    \[ RD = P(Y \mid X_{1},Z) - P(Y \mid X_{2},Z) \]

  • Relative Risk

\[ RR = \frac{P(Y \mid X_{1},Z)}{P(Y \mid X_{2},Z)} \]

  • Odds Ratio

\[ OR = \frac{P(Y \mid X_{1},Z)/(1-P(Y \mid X_{1},Z))}{P(Y \mid X_{2},Z)-(1-P(Y \mid X_{2},Z))} \]

Contoh Kasus:

Data menampilkan tabel kontingensi 3 arah, yang memperlihatkan hubungan pengaruh jenis kelamin, konsumsi minuman beralkohol terhadap perilaku seks bebas pada remaja. Data berjumlah 620 dengan proporsi 350 laki-laki, dan 270 perempuan.

\[\begin{array}{|c|c|c|c|} \hline \textbf{Jenis Kelamin} & \textbf{Alkohol} & \textbf{Seks Bebas: Ya} & \textbf{Seks Bebas: Tidak} \\ \hline Laki-Laki & Ya & 130 & 100 \\ Laki-Laki & Tidak & 30 & 90 \\ Perempuan & Ya & 70 & 90 \\ Perempuan & Tidak & 20 & 90 \\ \hline \end{array}\]

Tanpa memperthatikan variabel Z diperoleh tabel marginal sebagai berikut:

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Seks Bebas} & \textbf{Tidak Seks Bebas} & \textbf{Total} \\ \hline \textbf{Alkohol} & 200 & 190 & 390 \\ \textbf{Tidak Alkohol} & 50 & 180 & 230 \\ \hline \textbf{Total} & 250 & 370 & 620 \\ \hline \end{array}\]

Lakukan perhitungan manual:

  • Risk Difference:

    \[ RD = 0.51282 - 0.21739 = 0.29174 \]

  • Relative Risk:

    \[ RR = \frac{0.51282}{0.21739} = 2.35897 \]

  • Odds Ratio:

    \[ OR = \frac{200 \times 180}{190 \times 50} = 3.78947 \]

Dengan menggunakan syntax R:

# Data untuk laki-laki
a_laki <- 130  # Alkohol, Seks Bebas
b_laki <- 100  # Alkohol, Tidak Seks Bebas
c_laki <- 30   # Tidak Alkohol, Seks Bebas
d_laki <- 90   # Tidak Alkohol, Tidak Seks Bebas

# Menghitung proporsi risiko
P_ekspos_laki <- a_laki / (a_laki + b_laki)
P_tidak_ekspos_laki <- c_laki / (c_laki + d_laki)

# Menghitung ukuran asosiasi
RD_laki <- P_ekspos_laki - P_tidak_ekspos_laki
RR_laki <- P_ekspos_laki / P_tidak_ekspos_laki
OR_laki <- (a_laki * d_laki) / (b_laki * c_laki)

# Menampilkan hasil
cat("Beda Risiko (RD) Laki-laki: ", RD_laki, "\n")
## Beda Risiko (RD) Laki-laki:  0.3152174
cat("Risiko Relatif (RR) Laki-laki: ", RR_laki, "\n")
## Risiko Relatif (RR) Laki-laki:  2.26087
cat("Odds Ratio (OR) Laki-laki: ", OR_laki, "\n")
## Odds Ratio (OR) Laki-laki:  3.9

Ketiga ukuran asosiasi ini menggambarkan bahwa terdapat hubungan positif antara konsumsi alkohol dan perilaku seks bebas. Namun, perlu diperhatikan lagi bahwa hubungan ini tidak selalu bersifat kasual, dan faktor lain juga dapat berkontribusi terhadap perilaku seks bebas.

Jika, nilai Z diperhatikan, maka akan didapatkan ukuran asosiasi sebagai berikut:

  • Risk Difference:

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

  • Relative Risk

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

  • Odds Ratio

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

Misalkan kita akan menghitung ukuran asosiasi untuk orang mengonsumsi alkohol dibandingkan dengan yang tidak mengonsumsi alkohol terhadap perliaku seks bebas (untuk jenis kelamin laki-laki).

\[\begin{array}{|c|c|c|c|} \hline \textbf{Jenis Kelamin} & \textbf{Alkohol} & \textbf{Seks Bebas: Ya} & \textbf{Seks Bebas: Tidak} \\ \hline Laki-Laki & Ya & 130 & 100 \\ Laki-Laki & Tidak & 30 & 90 \\ Perempuan & Ya & 70 & 90 \\ Perempuan & Tidak & 20 & 90 \\ \hline \end{array}\]

Perhitungan manual:

  1. Kelompok Laki-Laki

    • Risk Difference :\[ RD = \frac{130}{130+100} - \frac{30}{30+90} = \frac{130}{230} - \frac{30}{120} = 0.3152174 \]

    • Relative Risk :

    \[ RR = \frac{\frac{130}{130+100}}{\frac{30}{30+90}} = 2.26087 \]

    • Odds Ratio :

    \[ OR = \frac{130 \times 90}{100 \times 30} = 3.9 \]

Dengan menggunakan syntax R:

# Data untuk laki-laki
a_laki <- 130  # Alkohol, Seks Bebas
b_laki <- 100  # Alkohol, Tidak Seks Bebas
c_laki <- 30   # Tidak Alkohol, Seks Bebas
d_laki <- 90   # Tidak Alkohol, Tidak Seks Bebas

# Menghitung proporsi risiko
P_ekspos_laki <- a_laki / (a_laki + b_laki)
P_tidak_ekspos_laki <- c_laki / (c_laki + d_laki)

# Menghitung ukuran asosiasi
RD_laki <- P_ekspos_laki - P_tidak_ekspos_laki
RR_laki <- P_ekspos_laki / P_tidak_ekspos_laki
OR_laki <- (a_laki * d_laki) / (b_laki * c_laki)

# Menampilkan hasil
cat("Beda Risiko (RD) Laki-laki: ", RD_laki, "\n")
## Beda Risiko (RD) Laki-laki:  0.3152174
cat("Risiko Relatif (RR) Laki-laki: ", RR_laki, "\n")
## Risiko Relatif (RR) Laki-laki:  2.26087
cat("Odds Ratio (OR) Laki-laki: ", OR_laki, "\n")
## Odds Ratio (OR) Laki-laki:  3.9

Interpretasi:

  • Laki-laki yang mengonsumsi alkohol memiliki peluang 0,3152174 lebih tinggi untuk melakukan seks bebas dibandingkan laki-laki yang tidak mengonsumsi alkohol.

  • Laki-laki yang mengonsumsi alkohol 2,26087 kali lebih mungkin melakukan seks bebas dibandingkan laki-laki yang tidak mengonsumsi alkohol.

  • Laki-laki yang mengonsumsi alkohol 3,9 kali lebih mungkin melakukan seks bebas dibandingkan laki-laki yang tidak mengonsumsi alkohol.

  1. Kelompok Perempuan

    • Risk Difference :\[ RD = \frac{70}{70+90} - \frac{20}{20+90} = \frac{70}{160} - \frac{20}{110} = 0.2556818 \]

    • Relative Risk :

      \[ RR = \frac{\frac{70}{70+90}}{\frac{20}{20+90}} = 2.40625 \]

    • Odds Ratio :

      \[ OR = \frac{70 \times 90}{20 \times 90} = 3.5 \]

Dengan menggunakan syntax R:

# Data untuk perempuan
a_perempuan <- 70   # Alkohol, Seks Bebas
b_perempuan <- 90   # Alkohol, Tidak Seks Bebas
c_perempuan <- 20   # Tidak Alkohol, Seks Bebas
d_perempuan <- 90   # Tidak Alkohol, Tidak Seks Bebas

# Menghitung proporsi risiko
P_ekspos_perempuan <- a_perempuan / (a_perempuan + b_perempuan)
P_tidak_ekspos_perempuan <- c_perempuan / (c_perempuan + d_perempuan)

# Menghitung ukuran asosiasi
RD_perempuan <- P_ekspos_perempuan - P_tidak_ekspos_perempuan
RR_perempuan <- P_ekspos_perempuan / P_tidak_ekspos_perempuan
OR_perempuan <- (a_perempuan * d_perempuan) / (b_perempuan * c_perempuan)

# Menampilkan hasil
cat("Beda Risiko (RD) Perempuan: ", RD_perempuan, "\n")
## Beda Risiko (RD) Perempuan:  0.2556818
cat("Risiko Relatif (RR) Perempuan: ", RR_perempuan, "\n")
## Risiko Relatif (RR) Perempuan:  2.40625
cat("Odds Ratio (OR) Perempuan: ", OR_perempuan, "\n")
## Odds Ratio (OR) Perempuan:  3.5

Interpretasi:

  • Perempuan yang mengonsumsi alkohol memiliki peluang 0,2556818 lebih tinggi untuk melakukan seks bebas dibandingkan perempuan yang tidak mengonsumsi alkohol.

  • Perempuan yang mengonsumsi alkohol 2,40625 kali lebih mungkin melakukan seks bebas dibandingkan perempuan yang tidak mengonsumsi alkohol.

  • Perempuan yang mengonsumsi alkohol 3,5 kali lebih mungkin melakukan seks bebas dibandingkan perempuan yang tidak mengonsumsi alkohol.

Conditional Independence

Terjadi ketika dua variabel menjadi independen setelah dikendalikan oleh variabel ketiga. Secara matematis, dua variabel X dan Y dikatakan independen secara kondisional terhadap variabel Z jika:

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

atau dalam bentuk frekuensi:

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

Salah satu untuk melihat kondisional independen adalah uji chi-square dengan hipotesis:

  • \(H_0\): Tidak ada hubungan antara konsumsi alkohol dan seks bebas setelah mengontrol jenis kelamin (kondisional independen).
  • \(H_A\): Ada hubungan antara konsumsi alkohol dan seks bebas bahkan setelah mengontrol jenis kelamin (tidak kondisional independen).

Dengan statistik uji:

\[ \chi^2 = \sum \frac{(O_{ij} - E_{ij})^2}{E_{ij}} \]

di mana: - \(O_{ij}\) = Frekuensi yang diamati dalam kategori (observed frequency). - \(E_{ij}\) = Frekuensi yang diharapkan dalam kategori (expected frequency), dihitung sebagai:

\[ E_{ij} = \frac{(\text{Total baris}) \times (\text{Total kolom})}{\text{Total keseluruhan}} \]

Derajat kebebasan (df) dihitung dengan:

\[ df = (r-1)(c-1) \]

di mana:

  • \(r\) = jumlah baris (kategori konsumsi alkohol).

  • \(c\) = jumlah kolom (kategori seks bebas).

Melakukan Perhitungan:

Tabel kontingensi untuk laki-laki:

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Seks Bebas (Ya)} & \textbf{Seks Bebas (Tidak)} & \textbf{Total} \\ \hline \textbf{Alkohol} & 130 & 100 & 230 \\ \textbf{Tidak Alkohol} & 30 & 90 & 120 \\ \hline \textbf{Total} & 160 & 190 & 350 \\ \hline \end{array}\]

Hitung Expected Frequency (\(E\)):

\[ E_{11} = \frac{(230 \times 160)}{350} = 105.14 \]

\[ E_{12} = \frac{(230 \times 190)}{350} = 124.86 \]

\[ E_{21} = \frac{(120 \times 160)}{350} = 54.86 \]

\[ E_{22} = \frac{(120 \times 190)}{350} = 65.14 \]

Hitung \(\chi^2\):

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

\[ \chi^2 = \frac{(130 - 105.14)^2}{105.14} + \frac{(100 - 124.86)^2}{124.86} + \frac{(30 - 54.86)^2}{54.86} + \frac{(90 - 65.14)^2}{65.14} \]

\[ \chi^2 = \frac{(24.86)^2}{105.14} + \frac{(-24.86)^2}{124.86} + \frac{(-24.86)^2}{54.86} + \frac{(24.86)^2}{65.14} \]

\[ \chi^2 = \frac{617.06}{105.14} + \frac{617.06}{124.86} + \frac{617.06}{54.86} + \frac{617.06}{65.14} \]

\[ \chi^2 = 5.87 + 4.94 + 11.25 + 9.48 \]

\[ \chi^2 = 31.54 \]

Tabel kontingensi untuk perempuan:

\[\begin{array}{|c|c|c|c|} \hline & \textbf{Seks Bebas (Ya)} & \textbf{Seks Bebas (Tidak)} & \textbf{Total} \\ \hline \textbf{Alkohol} & 70 & 90 & 160 \\ \textbf{Tidak Alkohol} & 20 & 90 & 110 \\ \hline \textbf{Total} & 90 & 180 & 270 \\ \hline \end{array}\]

Hitung Expected Frequency (\(E\)):

\[ E_{11} = \frac{(160 \times 90)}{270} = 53.33 \]

\[ E_{12} = \frac{(160 \times 180)}{270} = 106.67 \]

\[ E_{21} = \frac{(110 \times 90)}{270} = 36.67 \]

\[ E_{22} = \frac{(110 \times 180)}{270} = 73.33 \]

Hitung \(\chi^2\):

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

\[ \chi^2 = \frac{(70 - 53.33)^2}{53.33} + \frac{(90 - 106.67)^2}{106.67} + \frac{(20 - 36.67)^2}{36.67} + \frac{(90 - 73.33)^2}{73.33} \]

\[ \chi^2 = \frac{(16.67)^2}{53.33} + \frac{(-16.67)^2}{106.67} + \frac{(-16.67)^2}{36.67} + \frac{(16.67)^2}{73.33} \]

\[ \chi^2 = \frac{278.09}{53.33} + \frac{278.09}{106.67} + \frac{278.09}{36.67} + \frac{278.09}{73.33} \]

\[ \chi^2 = 5.21 + 2.61 + 7.58 + 3.79 \]

\[ \chi^2 = 19.19 \]

Dengan menggunakan syntax R:

data_laki <- matrix(c(130, 100, 30, 90), nrow = 2, byrow = TRUE,
                    dimnames = list(c("Alkohol_Yes", "Alkohol_No"), c("Seks_Bebas_Yes", "Seks_Bebas_No")))

data_perempuan <- matrix(c(70, 90, 20, 90), nrow = 2, byrow = TRUE,
                         dimnames = list(c("Alkohol_Yes", "Alkohol_No"), c("Seks_Bebas_Yes", "Seks_Bebas_No")))

# Uji Chi-Square untuk kelompok laki-laki
chi_laki <- chisq.test(data_laki, correct=FALSE)

# Menampilkan hasil uji Chi-Square
print("Hasil Uji Chi-Square untuk Laki-laki:")
## [1] "Hasil Uji Chi-Square untuk Laki-laki:"
print(chi_laki)
## 
##  Pearson's Chi-squared test
## 
## data:  data_laki
## X-squared = 31.574, df = 1, p-value = 1.92e-08
# Uji Chi-Square untuk kelompok perempuan
chi_perempuan <- chisq.test(data_perempuan, correct=FALSE)

# Menampilkan hasil uji Chi-Square
print("Hasil Uji Chi-Square untuk Perempuan:")
## [1] "Hasil Uji Chi-Square untuk Perempuan:"
print(chi_perempuan)
## 
##  Pearson's Chi-squared test
## 
## data:  data_perempuan
## X-squared = 19.176, df = 1, p-value = 1.192e-05

Interpretasi:

  • Tolak H0, maka dapat disimpulkan bahwa jenis kelamin (Z) mempengaruhi hubungan antara konsumsi alkohol (X) dan seks bebas (Y).

Inferensi Tabel Kontingensi Tiga Arah

Dalam analisis data kategorik, sering kali kita dihadapkan pada situasi di mana terdapat dua variabel utama yang ingin dianalisis hubungannya (misalnya variabel X dan Y), tetapi hubungan tersebut dapat dipengaruhi oleh variabel ketiga yang bersifat sebagai confounding atau variabel kontrol (Z). Untuk menangani situasi seperti ini, digunakanlah tabel kontingensi tiga arah, yang memungkinkan kita mengevaluasi asosiasi antara dua variabel sembari mengontrol efek dari variabel ketiga.

Sebagai contoh, kita ingin mengetahui apakah terdapat hubungan antara konsumsi alkohol (X) dengan kejadian perilau seks bebas (Y), dengan mempertimbangkan jenis kelamin sebagai variabel pengganggu (Z). Pendekatan yang digunakan adalah dengan membagi data ke dalam tabel parsial (yaitu tabel 2 × 2) untuk setiap strata dari variabel Z (misalnya, kelompok usia tertentu), serta membuat sebuah tabel marginal yang mengabaikan variabel Z.

Untuk menghitung kekuatan di tiap tabel parsial akan digunakan ukuran asosiasi yang umum digunakan, yakni Odds Ratio: \[ OR = \frac{a \cdot d}{b \cdot c} \] dengan:

  • \(a\) = jumlah kasus pada kategori positif Y dan positif X (misal: seks bebas dan alkohol)
  • \(b\) = jumlah kasus pada kategori positif Y dan negatif X (misal: seks bebas dan tidak alkohol)
  • \(c\) = jumlah kasus pada kategori negatif Y dan positif X (misal: tidak seks bebas dan alkohol)
  • \(d\) = jumlah kasus pada kategori negatif Y dan negatif X (misal: tidak seks bebas dan tidak alkohol)

Jika odds ratio parsial relatif konstan, kita dapat menghitung odds ratio bersama menggunakan estimasi Mantel-Haenszel.

\[ \text{OR}_{MH} = \frac{ \sum_{i=1}^{K} \frac{a_i d_i}{n_i} }{ \sum_{i=1}^{K} \frac{b_i c_i}{n_i} } \] Metode Mantel-Haenszel juga dapat digunakan jika ingin menguji apakah terdapat asosiasi antara variabel X dan Y dalam semua strata variabel Z.

Independensi Bersyarat Tabel Kontingensi 3 Arah

Dalam analisis data kategorik, Independensi Bersyarat merupakan bagian penting yang memungkinkan kita untuk memahami hubungan antara dua variabel utama dengan mempertimbangkan adanya variabel ketiga yang berpotensi mempengaruhi hubungan tersebut. Secara konseptual, independensi bersyarat terjadi ketika dua variabel X dan Y tidak saling bergantung satu sama lain pada setiap tingkat strata variabel ketiga (Z). Secara matematis dapat dituliskan sebagai berikut:

\[ OR(X,Y \mid Z) = 1 \]

Artinya Odds Ratio antara X dan Y adalah 1 pada setiap strata Z, yang berarti tidak adanya hubungan atau ketergantungan antara keduanya dalam konteks strata tersebut.

Sebagai contoh misalkan seorang peneliti ingin mengetahui apakah hubungan variabel X (konsumsi alkohol) dan Y (seks bebas) tetap signifikan dengan mempertimbangkan variabel Z (jenis kelamin) sebagai variabel kontrol. Jika OR =1 atau mendekati 1 pada setiap strata Z maka independensi bersyarat berlaku dan variabel X dan Y independensi bersyarat terhadap variabel Z (jenis kelamin), dan jika nilai OR bervariasi secara nyata pada setiap strata Z maka independensi bersyarat tidak berlaku.

Pengujian Statistik Independensi Besyarat dilakukan dengan uji Cochran-Mantel-Haenszel dengan hipotesis sebagai berikut:

Hipotesis:

  • \(H_0\) (Hipotesis nol): \(\theta_{XY(k)} = 1\) untuk semua \(k = 1, 2, \dots, K\)
  • \(H_1\) (Hipotesis alternatif): \(\theta_{XY(k)} \ne 1\) untuk setidaknya satu \(k\)

Statistik Uji CMH:

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

Keterangan:

  • \(n_{11k}\): Frekuensi aktual sel baris 1 kolom 1 pada strata ke-\(k\)
  • \(\mu_{11k}\): Ekspektasi frekuensi sel tersebut:

\[ \mu_{11k} = E(n_{11k}) = \frac{n_{1\cdot k} \cdot n_{\cdot 1k}}{n_{\cdot \cdot k}} \]

  • Varian dari \(n_{11k}\):

\[ \text{var}(n_{11k}) = \frac{n_{1\cdot k} \cdot n_{2\cdot k} \cdot n_{\cdot 1k} \cdot n_{\cdot 2k}}{n_{\cdot\cdot k}^2(n_{\cdot\cdot k}-1)} \]

Keputusan Uji:

Statistik uji CMH mengikuti distribusi chi-square dengan derajat kebebasan 1.

Tolak \(H_{0}\) jika,

\[ CMH > \chi^2_{(1)} \quad \text{atau p-value} < \alpha \]

Contoh Kasus:

\[\begin{array}{|c|c|c|c|} \hline \textbf{Jenis Kelamin} & \textbf{Alkohol} & \textbf{Seks Bebas: Ya} & \textbf{Seks Bebas: Tidak} \\ \hline Laki-Laki & Ya & 130 & 100 \\ Laki-Laki & Tidak & 30 & 90 \\ Perempuan & Ya & 70 & 90 \\ Perempuan & Tidak & 20 & 90 \\ \hline \end{array}\]

Hitung nilai ekspektasi:

\[ \mu_{11k} = \frac{n_{1\cdot k} \cdot n_{\cdot 1k}}{n_{\cdot\cdot k}} \]

\[\begin{align*} \mu_{11(1)} &= \frac{230 \times 160}{350} = 105.143 \\ \mu_{11(2)} &= \frac{160 \times 90}{270} = 53.33 \\ \end{align*}\]

Hitung nilai varians:

\[ \text{Var}(n_{11k}) = \frac{n_{1\cdot k} \cdot n_{2\cdot k} \cdot n_{\cdot 1k} \cdot n_{\cdot 2k}}{n_{\cdot\cdot k}^2 (n_{\cdot\cdot k} - 1)} \]

\[\begin{align*} \text{Var}(n_{11(1)}) &= \frac{230 \times 120 \times 160 \times 190}{350^2 \times 349} = 19.626 \\ \text{Var}(n_{11(2)}) &= \frac{160 \times 110 \times 90 \times 180}{270^2 \times 269} = 14.539 \\ \end{align*}\]

Hitung statistik uji CMH:

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

\[\begin{align*} X^2_{CMH} &= \frac{1724.227}{19.626 + 14.539} \\ &= \frac{1724.227}{34.165} \\ &= 50.468 \end{align*}\]

Nilai tersebut dibandingkan dengan distribusi \(\chi^2\) dengan derajat kebebasan 1.

Dengan \(\alpha = 0.05\), nilai kritis \(\chi^2_{1,0.05} = 3.841\). Karena \(174.087 > 3.841\), maka hipotesis nol ditolak, yang berarti ada hubungan signifikan antara konsumsi alkohol dan perilaku seks bebas setelah mengendalikan variabel jenis kelamin. Ini menunjukkan bahwa konsumsi alkohol dan perilaku seks bebas tidak bersifat conditional independence setelah jenis kelamin dikendalikan.

Dengan menggunakan syntax R:

# Membuat array 2x2x2 (2 lapisan untuk jenis kelamin)
data_cmh <- array(c(130, 30, 100, 90,  # Laki-laki
                    70, 20, 90, 90),   # Perempuan
                  dim = c(2, 2, 2), 
                  dimnames = list(
                    "Alkohol" = c("Ya", "Tidak"),
                    "Seks Bebas" = c("Ya", "Tidak"),
                    "Jenis Kelamin" = c("Laki-laki", "Perempuan")
                  ))

# Menampilkan tabel 2x2x2 untuk CMH
print("Tabel 2x2x2 untuk Uji CMH:")
## [1] "Tabel 2x2x2 untuk Uji CMH:"
print(data_cmh)
## , , Jenis Kelamin = Laki-laki
## 
##        Seks Bebas
## Alkohol  Ya Tidak
##   Ya    130   100
##   Tidak  30    90
## 
## , , Jenis Kelamin = Perempuan
## 
##        Seks Bebas
## Alkohol Ya Tidak
##   Ya    70    90
##   Tidak 20    90
# Melakukan Uji Cochran-Mantel-Haenszel (CMH)
library(vcd)
## Warning: package 'vcd' was built under R version 4.3.3
## Loading required package: grid
cmh_test <- mantelhaen.test(data_cmh, correct=FALSE)

# Menampilkan hasil uji CMH
print("Hasil Uji Cochran-Mantel-Haenszel:")
## [1] "Hasil Uji Cochran-Mantel-Haenszel:"
print(cmh_test)
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  data_cmh
## Mantel-Haenszel X-squared = 50.468, df = 1, p-value = 1.211e-12
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.565755 5.408009
## sample estimates:
## common odds ratio 
##             3.725

Didapatkan nilai p-value (1.211e-12) < 0.05 maka Tolak H0. Jadi, dapat disimpulkan bahwa terdapat hubungan signifikan antara konsumsi alkohol (X) dan seks bebas (Y) setelah mengontrol faktor jenis kelamin (Z).

Odds Ratio Bersama

Odds Ratio (OR) adalah ukuran asosiasi dalam tabel kontingensi 2x2 yang menunjukkan seberapa besar peluang suatu kejadian terjadi pada satu kelompok dibandingkan dengan kelompok lain. Namun, ketika data dikategorikan menjadi beberapa strata berdasarkan variabel perancu (seperti usia), kita tidak cukup hanya menghitung odds ratio dari satu tabel. Untuk itu, kita gunakan Odds Ratio Bersama atau Common Odds Ratio menggunakan pendekatan Mantel-Haenszel (MH).

\[ \hat{\theta}{MH} = \frac{\sum_{k=1}^{K} \left( \frac{n_{11k} n_{22k}}{n_{..k}} \right)}{\sum_{k=1}^{K} \left( \frac{n_{12k} n_{21k}}{n_{..k}} \right)} \]

Standar error log Odds Ratio bersama:

\[ \hat{\sigma}^{2}[\log(\hat{\theta}{MH})] = \frac{\sum (n{11k} + n_{12k})(n_{11k}n_{22k})/n_{..k}^{2}}{2(\sum n_{11k}n_{12k}/n_{..k})^{2}} + \frac{\sum(n_{11k} + n_{22k})(n_{11k} + n_{12k}) + (n_{12k} + n_{21k})(n_{11k} + n_{22k})]/n_{..k}^{2}}{2(\sum n_{11k}n_{12k}/n_{..k})^{2}} + \frac{\sum(n_{12k} + n_{21k})(n_{12k}n_{21k})/n_{..k}^{2}}{2(\sum n_{12k}n_{21k}/n_{..k})^{2}} \]

Interval Kepercayaan Odds Ratio Bersama:

\[ \log(\hat{\theta}{MH}) \pm Z{\alpha} \hat{\sigma}[\log(\hat{\theta}_{MH})] \]

Implementasi pada data konusmsi alkohol dan perilaku seks bebas:

\[ \begin{array}{|c|c|c|c|} \hline \textbf{Jenis Kelamin} & \textbf{Alkohol} & \textbf{Seks Bebas: Ya} & \textbf{Seks Bebas: Tidak} \\ \hline Laki-Laki & Ya & 130 & 100 \\ Laki-Laki & Tidak & 30 & 90 \\ Perempuan & Ya & 70 & 90 \\ Perempuan & Tidak & 20 & 90 \\ \hline \end{array} \]

Taksiran Odds Ratio Bersama:

\[\begin{align*} \text{Pembilang} &= \frac{(130)(90)}{350} + \frac{(70)(90)}{270} \\ &= 33.429 + 23.333 = 56.762 \\ \\ \text{Penyebut} &= \frac{(100)(30)}{350} + \frac{(90)(20)}{270} \\ &= 8.571 + 6.667 = 15.238 \\ \\ \hat{\theta}_{MH} &= \frac{56.762}{15.238} = 3.725 \end{align*}\]

Standar error log Odds Ratio:

Misalkan nilai standar error log Odds Ratio Bersama = 0.05.

Interval Kepercayaan Odds Ratio Bersama:

\[\begin{align*} \log(3.725) &\pm 1.96 \times 0.06 \\ 1.315 &\pm 0.098 \\ \Rightarrow \text{CI log} &= [1.217,\ 1.413] \\ \Rightarrow \text{CI OR} &= [e^{1.217},\ e^{1.413}] = [3.378,\ 4.108] \end{align*}\]

Interpretasi Hasil:

  • Odds ratio bersama menunjukkan bahwa konsumen alkohol 3,75 kali lebih mungkin melakukan seks bebas dibandingkan yang tidak mengonsumsi alkohol.

  • Standar error = 0.05 artinya variabilitas kecil pada estimasi cukup ratio

  • Interval Kepercayaan (3.378 hingga 4.108) menunjukkan bahwa kita cukup yakin bahwa odds ratio yang sebenarnya berada dalam interval ini.

  • Interval Kepercayaan tidak mencakup 1, hasil ini signifikan secara statistik, sekaligus mendukung adanya hubungan antara konsumsi alkohol dan perilaku seks bebas.

Dengan menggunakan syntax R:

# Membuat array 2x2x2 (2 lapisan untuk jenis kelamin)
library(epitools)
## 
## Attaching package: 'epitools'
## The following object is masked from 'package:vcd':
## 
##     oddsratio
data_cmh <- array(c(130, 30, 100, 90,  # Laki-laki
                    70, 20, 90, 90),   # Perempuan
                  dim = c(2, 2, 2), 
                  dimnames = list(
                    "Alkohol" = c("Ya", "Tidak"),
                    "Seks Bebas" = c("Ya", "Tidak"),
                    "Jenis Kelamin" = c("Laki-laki", "Perempuan")
                  ))

# Menampilkan tabel 2x2x2 untuk CMH
print("Tabel 2x2x2 untuk Uji CMH:")
## [1] "Tabel 2x2x2 untuk Uji CMH:"
print(data_cmh)
## , , Jenis Kelamin = Laki-laki
## 
##        Seks Bebas
## Alkohol  Ya Tidak
##   Ya    130   100
##   Tidak  30    90
## 
## , , Jenis Kelamin = Perempuan
## 
##        Seks Bebas
## Alkohol Ya Tidak
##   Ya    70    90
##   Tidak 20    90
# Melakukan Uji Cochran-Mantel-Haenszel (CMH)
cmh_test <- mantelhaen.test(data_cmh, correct=FALSE)

# Menampilkan hasil uji CMH
print("Hasil Uji Cochran-Mantel-Haenszel:")
## [1] "Hasil Uji Cochran-Mantel-Haenszel:"
print(cmh_test)
## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  data_cmh
## Mantel-Haenszel X-squared = 50.468, df = 1, p-value = 1.211e-12
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.565755 5.408009
## sample estimates:
## common odds ratio 
##             3.725

Uji Homogenitas Odds Ratio Statistik Breslow-Day

Dalam analisis tabel kontingensi 3 arah, kita ingin mengetahui apakah hubungan antara variabel X dan Y konsisten di seluruh strata dari variabel ketiga Z. Untuk mengetahui hal tersebut kita dapat menggunakan Uji Homogenitas Odds Ratio yang salah satu metode nya menggunakan statistik uji Breslow-Day.

Asosiasi dikatakan homogen apabila nilai odds ratio dari setiap strata memiliki nilai yang sama. Secara matematis kondisi ini dituliskan sebagai:

\[ \theta_{xy(1)} = \theta_{xy(2)} = \cdots = \theta_{xy(K)} \]

Jika kondisi ini terpenuhi maka dapat disimpulkan bahwa tidak terdapat interaksi atau hubungan X dan Y tidak bergantung pada variabel Z. Namun, jika kondisi tersebut tidak terpenuhi maka terindikasi bahwa hubungan X dan Y dipengaruhi oleh Z.

Untuk mengetahui hal tersebut akan dilakukan pengujian statistik menggunakan statistik uji Breslow-Day dengan hipotesis sebagai:

  • Hipotesis Nol (\(H_0\)):

    Tidak ada perbedaan odds ratio antar strata, yaitu:

    \[ \theta_{xy(1)} = \theta_{xy(2)} = \cdots = \theta_{xy(K)} \]

  • Hipotesis Alternatif (\(H_1\)):

    Setidaknya terdapat satu strata di mana nilai odds ratio berbeda dari yang lain.

Statistik Uji Breslow-Day:

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

Statistik uji Breslow-Day mengikuti distribusi chi square dengan derajat kebebasan = K-1

  • Jika nilai statistik uji lebih besar dari nilai kritis pada taraf signifikansi tertentu (misal 5%), maka kita menolak \(H_0\), yang berarti terdapat ketidakhomogenan odds ratio antar strata.
  • Sebaliknya, jika nilai statistik uji tidak signifikan, maka kita tidak menolak \(H_0\), sehingga asumsi homogenitas dapat diterima.

Implementasi pada contoh kasus hubungan konsumsi alkohol dan perilaku seks bebas:

\[ \begin{array}{|c|c|c|c|} \hline \textbf{Jenis Kelamin} & \textbf{Alkohol} & \textbf{Seks Bebas: Ya} & \textbf{Seks Bebas: Tidak} \\ \hline Laki-Laki & Ya & 130 & 100 \\ Laki-Laki & Tidak & 30 & 90 \\ Perempuan & Ya & 70 & 90 \\ Perempuan & Tidak & 20 & 90 \\ \hline \end{array} \]

Hitung Odds Ratio setiap strata:

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

Untuk setiap kelompok jenis kelamin:

\[\begin{align*} OR_1 &= \frac{(130 \times 90)}{(100 \times 30)} = \frac{11700}{3000} = 3.9 \\ OR_2 &= \frac{(70 \times 90)}{(90 \times 20)} = \frac{6300}{1800} = 3.5 \\ \end{align*}\]

Hitung Odds Ratio Bersama:

\[\begin{align*} \text{Pembilang} &= \frac{(130)(90)}{350} + \frac{(70)(90)}{270} \\ &= 33.429 + 23.333 = 56.762 \\ \\ \text{Penyebut} &= \frac{(100)(30)}{350} + \frac{(90)(20)}{270} \\ &= 8.571 + 6.667 = 15.238 \\ \\ \hat{\theta}_{MH} &= \frac{56.762}{15.238} = 3.725 \end{align*}\]

Hitung Ekspektasi \(\bar{a}_j\):

\[ -m_{j1} n_{j1} \hat{OR}{MH} + (n{j2} - m_{j1} + \hat{OR}{MH}(m{j1} + m_{j2})) \bar{a}j + (1 - \hat{OR}{MH}) \bar{a}_j^2 = 0 \]

\[\begin{align*} \bar{a}_1 &= \frac{m_{11} \cdot n_{11} \cdot \theta}{\theta(n_{11} + m_{21}) + m_{11} + n_{21}} \\ &= \frac{230 \cdot 160 \cdot 3.72}{3.72(160 + 120) + 230 + 190} \\ &= \frac{136896}{1461.6} \approx 93.7 \end{align*}\]

\[\begin{align*} \bar{a}_2 &= \frac{m_{1,2} \cdot n_{1,2} \cdot \theta}{\theta(n_{1,2} + m_{2,2}) + m_{1,2} + n_{2,2}} \\ &= \frac{160 \cdot 90 \cdot 3.72}{3.72(90 + 110) + 160 + 180} \\ &= \frac{53568}{744 + 340} = \frac{53568}{1084} \approx 49.42 \end{align*}\]

Hitung varians: \[ \text{Var}(a_j | \hat{OR}_{MH}) = \left( \frac{1}{a_j} + \frac{1}{b_j} + \frac{1}{c_j} + \frac{1}{d_j} \right)^{-1} \]

\[\begin{align*} \mathrm{Var}(a_1 \mid \hat{OR}_{MH}) &= \left( \frac{1}{a_1} + \frac{1}{b_1} + \frac{1}{c_1} + \frac{1}{d_1} \right)^{-1} \\ &= \left( \frac{1}{130} + \frac{1}{100} + \frac{1}{30} + \frac{1}{90} \right)^{-1} \\ &= \left( 0.00769 + 0.01 + 0.03333 + 0.01111 \right)^{-1} \\ &= \left( 0.06213 \right)^{-1} \approx 16.10 \end{align*}\]

\[\begin{align*} \mathrm{Var}(a_2 \mid \hat{OR}_{MH}) &= \left( \frac{1}{a_2} + \frac{1}{b_2} + \frac{1}{c_2} + \frac{1}{d_2} \right)^{-1} \\ &= \left( \frac{1}{70} + \frac{1}{90} + \frac{1}{20} + \frac{1}{90} \right)^{-1} \\ &= \left( 0.01429 + 0.01111 + 0.05 + 0.01111 \right)^{-1} \\ &= \left( 0.08651 \right)^{-1} \approx 11.56 \end{align*}\]

Hitung Statistik Uji Breslow-Day:

\[ \chi^2_{BD} = \sum_j \frac{(a_j - \bar{a}_j)^2}{\mathrm{Var}(a_j \mid \hat{OR}_{MH})} \]

\[\begin{align*} \chi^2_{BD} &= \frac{(130 - 93.7)^2}{16.10} + \frac{(70 - 56.3)^2}{11.56} \\&= \frac{(36.3)^2}{16.10} + \frac{(13.7)^2}{11.56} \\&= \frac{1317.69}{16.10} + \frac{187.69}{11.56} \\&\approx 81.86 + 16.23 = \boxed{98.09} \end{align*}\]

Terapkan Koreksi Tarone:

\[ X^2_{HBDT} = X^2_{HBD} - \frac{\sum_j (a_j - \bar{a}_j)^2}{\sum_j \text{Var}(a_j)} \]

\[ X^2_{HBDT} = 98.09 - \frac{1505.38}{27.66} = 98.09 \times 54.41 = 43.68 \]

Hitung P-Value:

\[ p = 1 - \chi^2(X^2_{HBDT}, df = 1) \]

\[ p = 1 - \chi^2(43.68, df = 1) \approx 3.87e-11 \]

Dengan menggunakan syntax R:

library(DescTools)
## Warning: package 'DescTools' was built under R version 4.3.3
data_bdt <- array(c(130, 30, 100, 90, 70, 20, 90, 90),
                  dim = c(2, 2, 2), 
                  dimnames = list(
                    "Alkohol" = c("Ya", "Tidak"),
                    "Seks Bebas" = c("Ya", "Tidak"),
                    "Jenis Kelamin" = c("Laki-laki", "Perempuan")
                  ))
print(data_bdt)
## , , Jenis Kelamin = Laki-laki
## 
##        Seks Bebas
## Alkohol  Ya Tidak
##   Ya    130   100
##   Tidak  30    90
## 
## , , Jenis Kelamin = Perempuan
## 
##        Seks Bebas
## Alkohol Ya Tidak
##   Ya    70    90
##   Tidak 20    90
breslow_test <- BreslowDayTest(data_bdt)
print(breslow_test)
## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  data_bdt
## X-squared = 0.078803, df = 1, p-value = 0.7789

Interpretasi:

Jika p-value < alpha (0.05), maka Tolak H0.

Jika p-value > alpha (0.05), maka Terima H0 artinya hubungan variabel X dan Y tidak dipengaruhi oleh variabel ketiga Z

Generelized Linear Model (GLM)

Generalized Linear Model (GLM) adalah sebuah kerangka statistik yang digunakan untuk memodelkan hubungan antara variabel dependen dengan satu atau lebih variabel independen. GLM adalah perluasan dari model linier klasik yang memungkinkan untuk menggunakan berbagai jenis distribusi untuk variabel dependen, bukan hanya distribusi normal.

GLM terdiri dari tiga komponen utama:

  1. Distribusi probabilitas untuk variabel dependen
    GLM memungkinkan penggunaan distribusi selain distribusi normal, seperti distribusi binomial, Poisson, dan lainnya. Distribusi ini tergantung pada sifat data yang dianalisis.

  2. Fungsi link
    Fungsi link menghubungkan nilai ekspektasi dari variabel dependen (\(E[Y]\)) dengan kombinasi linier dari variabel independen. Fungsi link ini mengatur bagaimana hubungan antara prediktor dan outcome.

  3. Kombinasi linier prediktor
    Kombinasi linier prediktor adalah bentuk linear dari variabel independen yang dihubungkan dengan koefisien model. Ini adalah bagian dari model yang memberikan prediksi dasar.

Exponential Family

Exponential family adalah sekumpulan distribusi probabilitas yang memiliki bentuk fungsi kepadatan atau fungsi massa probabilitas yang dapat dituliskan dalam bentuk eksponensial. Keluarga distribusi ini mencakup banyak distribusi penting dalam statistika, termasuk distribusi normal, binomial, Poisson, dan banyak lainnya.

Distribusi yang termasuk dalam exponential family memiliki bentuk umum sebagai berikut:

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

di mana:

  • \(\theta\) adalah natural parameter

  • \(\phi\) adalah dispersion parameter

  • \(a(\phi)\), \(b(\theta)\), dan \(c(y, \phi)\) adalah fungsi-fungsi tertentu yang berbeda untuk setiap distribusi

Distribusi yang termasuk dalam keluarga ini mencakup banyak distribusi umum seperti:

  • Distribusi Normal

  • Distribusi Binomial

  • Distribusi Poisson

  • Distribusi Gamma

Untuk menunjukkan bahwa suatu distribusi termasuk ke dalam exponential family, kita perlu menuliskan bentuk fungsi kepadatannya (PDF/PMF) dan menyusunnya agar sesuai dengan bentuk umum di atas.

Langkah-langkah:

  1. Tuliskan fungsi PDF atau PMF dari distribusi.
  2. Lakukan manipulasi aljabar agar dapat ditulis dalam bentuk: \[ f(y; \theta, \phi) = \exp\left( \frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi) \right) \]
  3. Identifikasi fungsi \(\theta\), \(a(\phi)\), \(b(\theta)\), dan \(c(y, \phi)\).

Contoh: Distribusi Poisson

Distribusi Poisson memiliki fungsi probabilitas:

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

Manipulasi ke Bentuk Exponential Family

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

Kita bisa tuliskan dalam bentuk:

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

Dengan substitusi \(\theta = \log \lambda\), maka \(\lambda = e^\theta\)

Identifikasi:

  • \(\theta = \log \lambda\)
  • \(b(\theta) = e^\theta\)
  • \(a(\phi) = 1\) (karena tidak ada parameter dispersi)
  • \(c(y, \phi) = -\log y!\)

Kesimpulan: Distribusi Poisson adalah anggota dari exponential family.

Regresi Logistik

Persamaan regresi logistik merupakan regresi linear, di mana nilai input dikombinasikan secara linear dengan sebuah bobot, yaitu menggunakan fungsi linear. Namun, regresi logistik membatasi hasil prediksi menjadi nilai antara 0 dan 1, dengan menggunakan fungsi aktivasi sigmoid.

Output yang diperoleh akan direpresentasikan dalam rentang nilai 0 hingga 1 dan sangatlah bersifat kurva S (S-shaped).

Regresi logistik menggambarkan hubungan antara satu atau lebih variabel independen dan menghasilkan nilai dalam dua kelas kategorik (biner). Model ini sangat cocok digunakan pada data dengan dua kelas keluaran (misalnya: ya/tidak, positif/negatif, sukses/gagal, dll).

Dengan kata lain, apabila di depan variabel bebas digunakan angka 1 mewakili kelas positif. Regresi logistik akan menghasilkan probabilitas untuk prediksi output kelas 1, dan variabel bebas tersebut akan dituliskan sebagai \(x = (x_1, x_2, ..., x_n)\).

Model secara umum berupaya memodelkan probabilitas bahwa \(Y = 1\) sebagai fungsi dari prediktor.

Beberapa contoh penerapan klasifikasi biner ini dapat ditemukan dalam berbagai kehidupan sehari-hari, seperti:

  • Memprediksi pelanggan akan tetap berlangganan dengan layanan produk digital, tetap menggunakan kartu kredit atau berhenti, dengan variabel-variabel seperti lama berlangganan, frekuensi penggunaan, status pembayaran, dan jumlah transaksi.
  • Memprediksi apakah email merupakan spam atau bukan dengan menggunakan fitur-fitur seperti kata-kata tertentu, penggunaan huruf kapital, dan struktur pesan.
  • Memprediksi apakah siswa akan lulus atau gagal ujian berdasarkan data akademik dan kehadiran mereka.
  • Menentukan apakah pasien terindikasi mengidap penyakit tertentu atau tidak berdasarkan hasil uji laboratorium dan gejala.
  • Analisis sentimen sosial: apakah suatu ulasan produk mengandung sentimen positif atau negatif.
  • Memprediksi apakah seseorang diterima di perguruan tinggi berdasarkan variabel prediktor seperti nilai rata-rata, nilai ujian masuk, dan partisipasi ekstrakurikuler.

Keunggulan Utama Regresi Logistik

Regresi logistik lebih mudah dibandingkan dengan pendekatan neural network karena modelnya lebih sederhana dan interpretasinya lebih jelas. Beberapa keuntungan lainnya adalah:

  • Tidak memerlukan asumsi distribusi normal pada variabel bebas.
  • Dapat digunakan untuk prediksi dan klasifikasi pada data biner.
  • Memungkinkan interpretasi langsung dari koefisien regresi sebagai odds ratio.
  • Cocok untuk data kategorik dan data yang tidak berdistribusi normal.
  • Hasil model mudah divisualisasikan dan dijelaskan.

Fungsi Aktivasi Sigmoid

Untuk mengubah output linear menjadi probabilitas, digunakan fungsi sigmoid. Fungsi sigmoid adalah:

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

Dengan fungsi ini, hasil prediksi akan bernilai antara 0 hingga 1. Beberapa interpretasi:

  • Jika hasil fungsi sigmoid lebih dari 0.5, maka output dianggap sebagai 1 (kelas positif).
  • Jika hasil fungsi kurang dari 0.5, maka output diklasifikasikan sebagai 0 (kelas negatif).
  • Jika mendekati 0.5, maka kemungkinan tidak terlalu pasti; bisa dilihat pada konteks atau dengan threshold berbeda.

Fungsi sigmoid inilah yang menjadikan regresi logistik berbeda dari model linear biasa.Plot kurva sigmoid sebagai berikut:

# Simulasi data untuk regresi logistik - Peluang Lulus Ujian
set.seed(123)
n <- 120
jam_belajar <- seq(0, 10, length.out = n)

# Ubah parameter regresi
intercept <- -3
slope <- 0.8

# Hitung log odds dan probabilitas
log_odds <- intercept + slope * jam_belajar
prob_lulus <- 1 / (1 + exp(-log_odds))

# Simulasi label (lulus atau tidak)
status_lulus <- rbinom(n, 1, prob_lulus)

# Buat data frame
data2 <- data.frame(
  jam_belajar = jam_belajar,
  status_lulus = status_lulus,
  prob_lulus = prob_lulus
)

# Visualisasi
plot(jam_belajar, status_lulus, pch = 16, col = "darkorange",
     xlab = "Jam Belajar", ylab = "Status Lulus (0 = Gagal, 1 = Lulus)",
     main = "Simulasi Regresi Logistik: Peluang Lulus Ujian")
lines(jam_belajar, prob_lulus, col = "darkgreen", lwd = 2)
abline(h = 0.5, col = "red", lty = 2)

legend("bottomright",
       legend = c("Data Biner", "Kurva Logistik", "Ambang 0.5"),
       col = c("darkorange", "darkgreen", "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 klasifkasi biner seperti deteksi penyakit, email spam, dan prediksi ya/tidak

Fungsi Link Logit dalam Regresi Logistik

Fungsi link (logit)

Fungsi link logit digunakan dalam regresi logistik untuk menghubungkan nilai rata-rata dari respon (\(\mu\)) dengan kombinasi linear dari prediktor. Fungsi ini didefinisikan sebagai:

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

Fungsi ini merepresentasikan log-odds dari probabilitas keberhasilan \(\mu\). Artinya, ia mengubah nilai probabilitas yang berada pada rentang \((0, 1)\) menjadi nilai yang tidak dibatasi dalam rentang \((-\infty, +\infty)\).

Model Regresi Logistik

Model regresi logistik secara umum dituliskan sebagai:

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

di mana: - \(\mu\) adalah probabilitas keberhasilan (misalnya, peluang sukses, lulus, positif, dll), - \(X\) adalah matriks variabel independen (fitur/prediktor), - \(\beta\) adalah vektor koefisien regresi.

Fungsi inverse link

Untuk memperoleh nilai probabilitas \(\mu\) dari kombinasi linear \(X\beta\), digunakan fungsi inverse dari logit, yaitu:

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

Fungsi ini dikenal sebagai fungsi sigmoid, yang memetakan nilai dari seluruh garis real ke interval \((0, 1)\), cocok untuk merepresentasikan probabilitas dalam regresi logistik.

Estimasi Parameter Dalam regresi logistik, tujuan utama dari proses estimasi adalah untuk menentukan parameter \(\beta\) yang memaksimalkan kesesuaian model terhadap data yang diamati. Karena model regresi logistik bukan linear dalam artian biasa, maka metode kuadrat terkecil (least squares) tidak dapat digunakan secara langsung.

Sebagai gantinya, digunakan metode Maximum Likelihood Estimation (MLE) atau Estimasi Likelihood Maksimum.

Prinsip dasar dari MLE adalah mencari nilai parameter \(\beta\) yang memaksimalkan fungsi likelihood, yaitu peluang bahwa model menghasilkan data seperti yang diamati.

Misalkan kita punya data pengamatan \((x_i, y_i)\), dengan: - \(y_i \in \{0, 1\}\) - \(\mu_i = P(y_i = 1 | x_i) = \frac{\exp(x_i^T \beta)}{1 + \exp(x_i^T \beta)}\)

Fungsi likelihood untuk seluruh data (dengan asumsi independen) dapat dituliskan sebagai:

\[ L(\beta) = \prod_{i=1}^{n} \mu_i^{y_i} (1 - \mu_i)^{1 - y_i} \]

Agar lebih mudah dihitung, digunakan log-likelihood:

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

Substitusi \(\mu_i = \frac{\exp(x_i^T \beta)}{1 + \exp(x_i^T \beta)}\), maka:

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

Proses estimasi dilakukan dengan mengoptimalkan fungsi log-likelihood ini, biasanya dengan metode numerik, seperti:

  • Newton-Raphson
  • Fisher Scoring
  • Iteratively Reweighted Least Squares (IRLS)
# 1. Simulasi data (ulang agar sesuai)
set.seed(123)
n <- 120
jam_belajar <- seq(0, 10, length.out = n)

# Parameter logit model
intercept <- -3
slope <- 0.8

# Hitung probabilitas lulus
log_odds <- intercept + slope * jam_belajar
prob_lulus <- 1 / (1 + exp(-log_odds))

# Simulasi label biner (0 = gagal, 1 = lulus)
status_lulus <- rbinom(n, 1, prob_lulus)

# Buat data frame
data_logistik <- data.frame(
  jam_belajar = jam_belajar,
  status_lulus = status_lulus
)

# 2. Estimasi model regresi logistik
model_logistik <- glm(status_lulus ~ jam_belajar, data = data_logistik, family = binomial)

# 3. Lihat hasil estimasi parameter
summary(model_logistik)
## 
## Call:
## glm(formula = status_lulus ~ jam_belajar, family = binomial, 
##     data = data_logistik)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.5339     0.5569  -4.550 5.35e-06 ***
## jam_belajar   0.7458     0.1326   5.626 1.85e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 155.387  on 119  degrees of freedom
## Residual deviance:  92.544  on 118  degrees of freedom
## AIC: 96.544
## 
## Number of Fisher Scoring iterations: 5
data_logistik$pred_prob <- predict(model_logistik, type = "response")
data_logistik$pred_logit <- predict(model_logistik, type = "link")

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
# --- 3. Plot Kurva Sigmoid (Probabilitas Lulus) ---
ggplot(data_logistik, aes(x = jam_belajar, y = status_lulus)) +
  geom_point(color = "gray60", alpha = 0.6, size = 2) +
  geom_line(aes(y = pred_prob), color = "blue", size = 1.2) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "red") +
  labs(title = "Kurva Sigmoid Regresi Logistik",
       x = "Jam Belajar",
       y = "Probabilitas / Status Lulus") +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Evaluasi Model

# --- 1. Prediksi klasifikasi dengan threshold 0.5 ---
data_logistik$pred_class <- ifelse(data_logistik$pred_prob > 0.5, 1, 0)

# Confusion Matrix
table_pred <- table(Prediksi = data_logistik$pred_class, Aktual = data_logistik$status_lulus)
print(table_pred)
##         Aktual
## Prediksi  0  1
##        0 31 10
##        1 11 68
# Hitung akurasi
akurasi <- sum(diag(table_pred)) / sum(table_pred)
cat("Akurasi model:", round(akurasi, 4), "\n")
## Akurasi model: 0.825
# --- 2. Nilai AIC dan Deviance ---
cat("AIC:", AIC(model_logistik), "\n")
## AIC: 96.5443
cat("Null Deviance:", model_logistik$null.deviance, "\n")
## Null Deviance: 155.3872
cat("Residual Deviance:", model_logistik$deviance, "\n")
## Residual Deviance: 92.5443
# --- 3. Pseudo R-squared (McFadden's R²) ---
pseudoR2 <- 1 - (model_logistik$deviance / model_logistik$null.deviance)
cat("Pseudo R² (McFadden):", round(pseudoR2, 4), "\n")
## Pseudo R² (McFadden): 0.4044

Regresi Poisson

Regresi Poisson adalah salah satu model dalam keluarga Generalized Linear Model (GLM) yang digunakan untuk memodelkan data diskrit berupa hitungan (jumlah kejadian dalam suatu interval waktu, ruang, atau unit lain).

Contoh kasus:

  • Jumlah pelanggan yang datang ke toko dalam 1 hari

  • Jumlah kecelakaan di suatu jalan dalam sebulan

  • Jumlah email spam yang diterima per minggu

Fungsi link: \[ \log(\lambda) = \beta_0 + \beta_1x \]

Simulasi Data dan Evaluasi Model:

set.seed(1)
n <- 100
x <- rnorm(n, mean = 0, sd = 1)
lambda <- exp(0.5 + 0.4 * x)
y <- rpois(n, lambda)
data_poisson <- data.frame(x = x, y = y)

model_poisson <- glm(y ~ x, data = data_poisson, family = poisson)
summary(model_poisson)
## 
## Call:
## glm(formula = y ~ x, family = poisson, data = data_poisson)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.30986    0.09134   3.392 0.000693 ***
## x            0.50125    0.09009   5.564 2.64e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 137.41  on 99  degrees of freedom
## Residual deviance: 106.15  on 98  degrees of freedom
## AIC: 303.81
## 
## Number of Fisher Scoring iterations: 5
exp(coef(model_poisson)) 
## (Intercept)           x 
##    1.363235    1.650778
model_poisson$deviance / df.residual(model_poisson)
## [1] 1.083182

Jika nilai deviance/df > 1.5, kemungkinan terjadi overdispersion.

Visualisasi

data_poisson$pred <- predict(model_poisson, type = "response")

plot(data_poisson$x, data_poisson$y, pch = 16, col = "gray",
     main = "Data dan Hasil Prediksi",
     xlab = "x", ylab = "y")

ord <- order(data_poisson$x)
lines(data_poisson$x[ord], data_poisson$pred[ord], col = "blue", lwd = 2)

Inferensi GLM

Mencari Ekspektasi dan Variansi dalam GLM

Ekspektasi

Jika diturunkan berdasarkan fungsi momen:

\[ \mathbb{E}(Y) = \int y f(y; \theta) dy = \mu \]

Untuk keluarga eksponensial:

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

atau:

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

Maka ekspektasi turunan pertama:

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

Dan ekspektasi turunan pertama:

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

Maka:

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

Varians

Turunan kedua:

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

Sehingga:

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

Metode Penaksiran Parameter

Maximum Likelihood Estimation (MLE)

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

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

Metode Optimisasi

Newton-Raphson

  • Menggunakan score vector (gradien)
  • Menggunakan Hessian matrix

Iterasi:

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

Fisher Scoring

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

IRLS (Iteratively Reweighted Least Square)

  • Modifikasi dari Fisher scoring, hasil estimasi mirip dengan Least Square.

Implementasi Newton-Raphson

Statistik score ke-j:

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

Turunan kedua:

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

Taylor expansion:

\[ U(\beta') \approx U(\beta) + H(\beta)(\beta' - \beta) \]

Estimasi parameter:

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

Diagnostik Model GLM

Diagnostik digunakan untuk mengevaluasi apakah model sudah tepat.

  • Uji formal
  • Grafik antara nilai prediksi vs nilai aktual

Statistik Devian

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

Devian adalah:

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

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

Statistik Chi-Kuadrat Pearson

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

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

  • Jika signifikan → model lebih baik daripada tanpa model.

Detai 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_i(\boldsymbol{x}) = \frac{\exp(\beta_0 + \beta_1 x_i)}{1 + \exp(\beta_0 + \beta_1 x_i)} \]

Log-likelihood untuk n observasi:

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

Estimasi dengan Newton-Raphson

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

Fungsi Log-Likelihood

Model regresi logistik untuk probabilitas:

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

Log-likelihood untuk n observasi:

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

Langkah-Langkah Newton-Raphson

  1. Turunan Pertama (Score Function):

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

  1. Turunan Kedua (Hessian Matrix):

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

  1. Iterasi Newton-Raphson:

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

Estimasi MLE

# Data simulasi
set.seed(1)
x <- rnorm(100)
beta_true <- c(-1.5, 2.2)  
X <- cbind(1, x)
eta <- X %*% beta_true
p <- 1 / (1 + exp(-eta))
y <- rbinom(100, 1, p)

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

for (i in 1:max_iter) {
    eta <- X %*% beta
    pi_hat <- 1 / (1 + exp(-eta))
    W <- diag(as.numeric(pi_hat * (1 - pi_hat)))
    z <- eta + solve(W) %*% (y - pi_hat)
    beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
    if (sum(abs(beta_new - beta)) < 0.0001) break 
    beta <- beta_new
}
beta
##        [,1]
##   -2.295432
## x  2.561043

Uji Wald

# Data simulasi
set.seed(123)
x <- rnorm(150)  # Sample size modified
beta_true <- c(-0.8, 1.5)  # Parameters modified
X <- cbind(1, x)
eta <- X %*% beta_true
p <- 1 / (1 + exp(-eta))
y <- rbinom(150, 1, p)

# Fit model
model <- glm(y ~ x, family = binomial)
beta_hat <- coef(model)["x"]
se_beta <- summary(model)$coefficients["x", "Std. Error"]  # Fixed space

Z <- beta_hat / se_beta
Wald_stat <- Z^2
p_value <- 1 - pchisq(Wald_stat, df = 1)

cat("Wald Statistic:", round(Wald_stat, 3), 
    "\np-value:", format.pval(p_value, digits = 3))
## Wald Statistic: 27.873 
## p-value: 1.3e-07

Jika p-value < alpha (0.05), maka Terima Hipotesis nol.

Detail Metode Estimasi dan Inferensi Regresi Poisson

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

Model Regresi Poisson

Fungsi distribusi Poisson:

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

Model regresi Poisson:

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

Estimasi Parameter (MLE)

Fungsi log-likelihood:

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

Dengan:

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

Estimasi menggunakan metode iterasi (IRLS)

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

Tahap 1: Definisikan model regresi Poisson:

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

sehingga:

\[ \lambda_i = \exp(\eta_i) \]

Tahap 2: Mencari log-likelihood yang dimaksimalkan:

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

Tahap 3: Formulasi iteratif IRLS:

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

Dengan:

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

dan

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

fEstimasi IRLS Manual

# Fungsi untuk estimasi manual regresi Poisson
poisson_irls <- function(X, y, max_iter = 50, tol = 1e-6) {
  beta <- rep(0, ncol(X))  # Inisialisasi parameter
  
  for (iter in 1:max_iter) {
    eta <- X %*% beta
    lambda <- exp(eta)
    W <- diag(as.vector(lambda))  # Matriks bobot
    z <- eta + (y - lambda) / lambda  # Variabel respon adjusted
    
    beta_new <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% z
    
    if (max(abs(beta_new - beta)) < tol) break
    beta <- beta_new
  }
  
  list(coefficients = beta, iterations = iter)
}

##Contoh
set.seed(123)
n <- 200
X <- cbind(1, rnorm(n))  # Matriks desain dengan intercept
beta_true <- c(0.5, 0.8)
y <- rpois(n, exp(X %*% beta_true))

# Estimasi parameter
result <- poisson_irls(X, y)
cat("Estimasi parameter:\n")
## Estimasi parameter:
print(result$coefficients)
##           [,1]
## [1,] 0.5095712
## [2,] 0.7923024
cat("Jumlah iterasi:", result$iterations)
## Jumlah iterasi: 9

Inferensi Statistik Uji Wald

# Menggunakan model dari glm() untuk contoh
model <- glm(y ~ X[,-1], family = poisson)

# Statistik Wald
wald_test <- function(model) {
  summ <- summary(model)
  coefs <- summ$coefficients
  z <- coefs[, "Estimate"] / coefs[, "Std. Error"]
  p_value <- 2 * pnorm(-abs(z))
  
  data.frame(
    Estimate = coefs[, "Estimate"],
    Std.Error = coefs[, "Std. Error"],
    z.value = z,
    p.value = p_value
  )
}

wald_results <- wald_test(model)
print(wald_results)
##              Estimate  Std.Error   z.value      p.value
## (Intercept) 0.5095711 0.05885608  8.657918 4.804596e-18
## X[, -1]     0.7923019 0.04433252 17.871800 1.955815e-71

Evaluasi Model

# AIC dan BIC
cat("AIC:", AIC(model), "\n")
## AIC: 635.3746
cat("BIC:", BIC(model), "\n")
## BIC: 641.9713
cat("Devians:", deviance(model))
## Devians: 216.0077

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:

Simulasi Data

# Simulasi variabel
n <- 500

set.seed(123)
job_type <- sample(c("Freelance", "Employee"), n, replace = TRUE)  # nominal
satisfaction <- sample(c("Very Unsatisfied", "Unsatisfied", "Neutral", "Satisfied", "Very Satisfied"),
                       n, replace = TRUE, prob = c(0.1, 0.2, 0.3, 0.25, 0.15))  # ordinal
work_hours <- rnorm(n, mean = 40, sd = 10)  # rasio

# Logit model: asumsi dummy untuk Red, urutan satisfaction sebagai ordinal, dan jam kerja numerik
logit_p <- -1.5 +
           0.6 * (job_type == "Freelance") +
           0.7 * as.numeric(factor(satisfaction, 
                                   levels = c("Very Unsatisfied", "Unsatisfied", "Neutral", "Satisfied", "Very Satisfied"), 
                                   ordered = TRUE)) +
           0.04 * work_hours

p <- 1 / (1 + exp(-logit_p))
set.seed(57)
success <- rbinom(n, 1, p)

library(tibble)
## Warning: package 'tibble' was built under R version 4.3.3
sim_data <- tibble(success, job_type, satisfaction, work_hours)

head(sim_data)

Interpretasi: Dataset berisikan 500 observasi dengan variabel warna mobil, kepuasan, jam kerja, dan status keberhasilan.

Eksplorasi Data

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 %>%
  group_by(success) %>%
  summarise(
    Jumlah = n(),
    Rata2_JamKerja = mean(work_hours),
    .groups = "drop"
  )

Interpretasi: Distribusi jumlah sukses dan tidak sukses adalah 445 dan 55. Rata-rata jam kerja pada masing-masing grup adalah 35.28656 dan 40.55633.

Perlakuan Variabel Ordinal

Treat Sebagai Nominal (Dummy)

sim_data_nominal <- sim_data %>%
mutate(
satisfaction = factor(satisfaction, levels = c("Very Unsatisfied", "Unsatisfied", "Neutral", "Satisfied", "Very Satisfied"))
)
model_nominal <- glm(success ~ job_type + satisfaction + work_hours, data = sim_data_nominal, family = binomial)
summary(model_nominal)
## 
## Call:
## glm(formula = success ~ job_type + satisfaction + work_hours, 
##     family = binomial, data = sim_data_nominal)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -1.78678    0.73609  -2.427 0.015207 *  
## job_typeFreelance           0.83875    0.31746   2.642 0.008240 ** 
## satisfactionUnsatisfied     0.56447    0.45251   1.247 0.212240    
## satisfactionNeutral         1.25094    0.46212   2.707 0.006790 ** 
## satisfactionSatisfied       1.90797    0.52677   3.622 0.000292 ***
## satisfactionVery Satisfied  3.63136    1.07618   3.374 0.000740 ***
## work_hours                  0.06026    0.01584   3.804 0.000143 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 346.52  on 499  degrees of freedom
## Residual deviance: 294.61  on 493  degrees of freedom
## AIC: 308.61
## 
## Number of Fisher Scoring iterations: 7

Keterangan: Model ini menggunakan satisfaction sebagai variabel dummy, baseline ““Very Unsatisfied”. Setiap koefisien membandingkan kategori terhadap baseline.

Interpretasi Koefisisen

Intercept (-1.78678)

  • Ini adalah log-odds dasar untuk individu yang memiliki warna mobil selain Blue (misalnya Black), tingkat kepuasan Very Unsatisfied, dan jumlah jam kerja = 0.
  • p = 0.0152 (signifikan pada level 5%), artinya intercept berbeda signifikan dari 0.
  • Odds Ratio = exp(-1.78678) = 0.17 → Peluang sukses hanya sekitar 17% dari baseline (sangat rendah).

car_colorBlue (0.83875)

  • Individu dengan warna mobil Blue memiliki log-odds sukses lebih tinggi sebesar 0.839 dibanding warna referensi (misalnya Black).
  • p = 0.0082 (signifikan), menunjukkan bahwa perbedaan ini bermakna secara statistik.
  • Odds Ratio = exp(0.83875) = 2.31 → Peluang sukses pengguna mobil Blue sekitar 2.3 kali lebih tinggi dibanding pengguna mobil Black.

satisfactionUnsatisfied (0.56447)

  • Individu dengan kepuasan Unsatisfied memiliki log-odds lebih tinggi sebesar 0.564 dibanding yang Very Unsatisfied.
  • p = 0.2122 (tidak signifikan), sehingga peningkatan ini tidak bisa dianggap bermakna.
  • Odds Ratio = exp(0.56447) = 1.76 → Peluang sukses meningkat 1.76 kali, namun tidak signifikan.

satisfactionNeutral (1.25094)

  • Individu dengan kepuasan Netral memiliki log-odds lebih tinggi sebesar 1.25 dibanding yang Very Unsatisfied.
  • p = 0.0068 (signifikan pada level 1%), menunjukkan peningkatan yang bermakna.
  • Odds Ratio = exp(1.25094) = 3.49 → Peluang sukses meningkat sekitar 3.5 kali dibanding Very Unsatisfied.

satisfactionSatisfied (1.90797)

  • Individu dengan kepuasan Satisfied memiliki log-odds lebih tinggi sebesar 1.91 dibanding Very Unsatisfied.
  • p = 0.0003 (sangat signifikan).
  • Odds Ratio = exp(1.90797) = 6.74 → Peluang sukses hampir 7 kali lipat dibanding individu yang sangat tidak puas.

work_hours (0.06026)

  • Setiap tambahan 1 jam kerja per minggu meningkatkan log-odds sebesar 0.06026.
  • p = 0.0001 (sangat signifikan).
  • Odds Ratio = exp(0.06026) = 1.062 → Setiap tambahan 1 jam kerja meningkatkan peluang sukses sebesar 6.2%.

Interpretasi Goodness-of-Fit

  • Null Deviance (346.52): Deviance model tanpa prediktor

  • Residual deviance (294.61): Deviance model dengan prediktor

    • Penurunan nilai deviance dari null deviance ke residual deviance menunjukkan model membawa informasi yang cukup baik.
  • AIC (308.61): Digunakan untuk membandingkan model, semakin kecil AIC maka semakin baik model dalam menyeimbangkan goodness-of-fit dan kompleksitas.

Signifikansi Model

  • Variabel satisfaction secara keseluruhan signifikan dalam meningkatkan peluang sukses. Terutama kategori Neutral, Satisfied, dan Very Satisfied, yang memiliki pengaruh signifikan (p < 0.01).

  • Warna mobil Blue juga menunjukkan pengaruh signifikan terhadap peluang sukses (p = 0.0082), dengan peluang sukses sekitar 2.3 kali lebih tinggi dibanding warna referensi.

  • Jumlah jam kerja (work_hours) signifikan secara statistik (p < 0.001), artinya semakin banyak jam kerja, semakin tinggi peluang sukses individu.

  • Kategori satisfaction Unsatisfied memiliki efek positif namun tidak signifikan (p = 0.2122), sehingga tidak cukup bukti bahwa kategori ini berbeda secara signifikan dari Very Unsatisfied.

Kesimpulan Praktis

  • Tingkat kepuasan individu merupakan prediktor paling kuat terhadap keberhasilan. Semakin tinggi tingkat kepuasan, semakin besar peluang sukses, dengan peningkatan hingga 38 kali lipat untuk kategori Very Satisfied dibanding Very Unsatisfied.

  • Warna mobil biru secara mengejutkan memiliki hubungan positif terhadap keberhasilan, namun ini bisa jadi artefak dari data simulasi. Interpretasi ini perlu dikaji lebih lanjut dengan variabel pengganti yang lebih logis.

  • Jam kerja memberikan kontribusi positif yang signifikan terhadap keberhasilan, dengan tambahan 1 jam kerja meningkatkan peluang sukses sebesar 6.2%.

  • Meskipun model ini sudah menjelaskan variabilitas data lebih baik daripada model null (lihat deviance dan AIC), masih ada ruang untuk perbaikan, misalnya dengan menambahkan variabel yang lebih bermakna seperti jenis pekerjaan atau motivasi pribadi.

Treat Sebagai Rasio (Numeric Rank)

library(dplyr)
sim_data_numeric <- sim_data %>%
  mutate(
    satisfaction_numeric = case_when(
      satisfaction == "Very Unsatisfied" ~ 1,
      satisfaction == "Unsatisfied" ~ 2,
      satisfaction == "Neutral"  ~ 3,
      satisfaction == "Satisfied" ~ 4,
      satisfaction == "Very Satisfied" ~ 5
  )
)
model_numeric <- glm(success ~ job_type + satisfaction_numeric + work_hours, data = sim_data_numeric, family = binomial)
summary(model_numeric)
## 
## Call:
## glm(formula = success ~ job_type + satisfaction_numeric + work_hours, 
##     family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -2.68665    0.75758  -3.546 0.000391 ***
## job_typeFreelance     0.86532    0.31581   2.740 0.006144 ** 
## satisfaction_numeric  0.72901    0.14034   5.195 2.05e-07 ***
## work_hours            0.06099    0.01590   3.836 0.000125 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 346.52  on 499  degrees of freedom
## Residual deviance: 295.91  on 496  degrees of freedom
## AIC: 303.91
## 
## Number of Fisher Scoring iterations: 6

Interpretasi: Model menggunakan satisfaction sebagai angka bertingkat 1-5. Koefisien mengukur kenaikan satu tingkat kepuasan terhadap peluang sukses.

Perbandingan Model

list(
AIC_Nominal = AIC(model_nominal),
AIC_Numeric = AIC(model_numeric)
)
## $AIC_Nominal
## [1] 308.6077
## 
## $AIC_Numeric
## [1] 303.9145

Keterangan: Model dengan nilai AIC yang lebih kecil cenderung lebih disukai.

Interpretasi Koefisien Intercept (-2.68665) Ini adalah log-odds dasar untuk individu dengan tipe pekerjaan Employee, kepuasan kerja terendah (Very Unsatisfied = 1), dan jam kerja = 0. Tidak dapat diinterpretasikan secara praktis karena tidak realistis (jam kerja = 0), namun tetap menjadi baseline referensi.

job_typeFreelance (0.86532) Individu dengan tipe pekerjaan Freelance memiliki log-odds peluang sukses yang lebih tinggi sebesar 0.865 dibanding Employee.

  • p-value = 0.006144 (signifikan)

  • Odds Ratio = exp(0.86532) = 2.376 → Peluang sukses freelancer sekitar 2.4 kali lebih tinggi dibandingkan pegawai.

satisfaction_numeric (0.72901) Setiap peningkatan satu tingkat kepuasan kerja (misalnya dari Unsatisfied ke Neutral) meningkatkan log-odds sukses sebesar 0.729.

  • p-value < 0.001 (sangat signifikan)

  • Odds Ratio = exp(0.72901) = 2.073 → Setiap kenaikan tingkat kepuasan kerja meningkatkan peluang sukses sebesar 107%.

work_hours (0.06099) Setiap penambahan satu jam kerja per minggu meningkatkan log-odds sukses sebesar 0.061.

  • p-value = 0.000125 (signifikan)

  • Odds Ratio = exp(0.06099) = 1.063 → Setiap tambahan satu jam kerja per minggu meningkatkan peluang sukses sekitar 6.3%.

Interpretasi Goodness-of-Fit - Null deviance (346.52): Deviance model tanpa prediktor.

  • Residual deviance (295.91): Deviance model dengan prediktor.→ Penurunan deviance menunjukkan bahwa prediktor memberikan tambahan informasi yang signifikan.

  • AIC (303.91): Digunakan untuk membandingkan model — semakin kecil, semakin baik.

Signifikansi Model Variabel satisfaction_numeric, work_hours, dan job_typeFreelance semuanya signifikan dalam meningkatkan peluang sukses.

Model secara umum baik dalam menjelaskan variabilitas respons (success).

Kesimpulan - Kepuasan kerja dan jam kerja memiliki pengaruh positif terhadap peluang sukses.

  • Pekerja freelance memiliki peluang lebih tinggi dibandingkan pekerja tetap.

  • Model yang dibangun secara statistik signifikan dan mampu menjelaskan variansi data.

Goodeness-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.1497989 (df=7)
## 
## $McFadden_R2_Numeric
## 'log Lik.' 0.1460278 (df=4)

Interpretasi: McFadden’s R^2 mengukur goodness-of-fit. Semakin besar nilainya, semakin baik model memprediksi dibandingkan model null.

Visualisasi Prediksi

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

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 = work_hours, y = predicted, color = as.factor(satisfaction_numeric))) +
geom_point(alpha = 0.6) +
labs(title = "Prediksi Probabilitas (Ordinal sebagai Numeric)", x = "Work Hours", y = "Prediksi Probabiliti")

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

Interpretasi: Visualisasi hubungan antara pendapatan dan probabilitas sukses berdasarkan tingkat kepuasan dengan dua pendekatan ordinal. Secara numerik, memperhatikan urutan dari Very Unsatisfied hingga Very Satisfied. Namun, secara nominal tidak memperhatikan urutan tingkat kepuasan nya.

Ringkasan Koefisien Model

library(broom)
## Warning: package 'broom' was built under R version 4.3.3
library(knitr)
## Warning: package 'knitr' was built under R version 4.3.3
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.3.3
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
# Ringkasan model nominal
tidy(model_nominal) %>%
  kable(format = "latex", booktabs = TRUE,
        caption = "Ringkasan Koefisien Model dengan Ordinal sebagai Nominal") %>%
  kable_styling(latex_options = c("hold_position","striped"))
# Ringkasan model numeric
tidy(model_numeric) %>%
kable(format = "latex", booktabs = TRUE, caption = "Ringkasan Koefisien Model dengan Ordinal sebagai Numeric") %>% kable_styling(latex_options = c("hold_position", "striped"))

Interpretasi: Tabel memperlihatkan estimasi koefisien, standar error, nilai z, dan p-value. Koefisien positif meningkatkan log-odds keberhasilan, koefisien negatif akan menurunkan log-odds keberhasilan.

Kesimpulan

  • Job Type: Freelance memiliki peluang sukses lebih tinggi dibandingkan dengan employee

  • Satisfaction:

    • Jika dianggap dummy (nominal) : tiap tingkat akan dibandingkan langsung dengan tingkat Very Unsatisfied

    • Jika dianggap numeric : tiap kenaikan satu tingkat kepuasan meningkatkan peluang sukses.

  • Work Hours: Waktu bekerja yang lebih tinggi meningkatkan peluang sukses.

Pada akhirnya, model dengan perlakuan ordinal sebagai numerik akan memiliki interpretasi lebih sederhana jika asumsi jarak antar tingkatannya konsisten.

Pemilihan Model Regresi Logistik dan Evaluasi

Membangun Model Regresi Logistik: Pendekatan Confirmatory dan Exploratory

Dalam analisis regresi logistik, pemilihan model yang tepat sangat penting untuk memperoleh model yang dapat memprediksi probabilitas kejadian suatu peristiwa (respon biner) secara akurat. Dua pendekatan utama yang sering digunakan dalam membangun model adalah pendekatan Konfirmatori dan Eksploratori.

1. Konfirmatori (Pendekatan Konfirmatori)

Pendekatan ini digunakan apabila peneliti sudah memiliki teori atau hipotesis yang jelas mengenai hubungan antara variabel prediktor dengan variabel respon.

Ciri-ciri:

  • Model dikembangkan berdasarkan teori atau temuan dari penelitian sebelumnya.
  • Fokus utama adalah menguji apakah efek dari variabel signifikan secara statistik, bukan hanya sekadar mencari model terbaik.
  • Peneliti biasanya telah menentukan model terlebih dahulu, lalu mengujinya dengan data, termasuk kemungkinan adanya interaksi atau efek tambahan lainnya.
  • Pengujian dilakukan dengan metode statistik seperti Uji Likelihood Ratio Test untuk menilai apakah efek dari variabel signifikan tanpa adanya efek lain.

Contoh penggunaan:

Misalnya, berdasarkan teori diasumsikan bahwa faktor x1 dan x2 mempengaruhi kemungkinan seseorang membeli produk. Maka, model regresi logistik akan langsung mengikutkan x1 dan x2 untuk diuji apakah kontribusinya signifikan.

2. Eksploratori (Pendekatan Eksploratori)

Pendekatan ini digunakan ketika peneliti belum memiliki teori yang jelas, dan bertujuan untuk mengeksplorasi hubungan antara variabel.

Ciri-ciri:

  • Model dikembangkan secara bertahap dengan tujuan menemukan kombinasi prediktor terbaik.
  • Pemilihan variabel dilakukan dengan kriteria statistik seperti AIC, deviance, atau log-likelihood.

Proses pemilihan model:

  • Forward Selection: Dimulai dari model kosong, kemudian variabel ditambahkan satu per satu jika signifikan.
  • Backward Elimination: Dimulai dari model penuh, kemudian variabel yang tidak signifikan dihapus.
  • Stepwise Selection: Kombinasi antara penambahan dan penghapusan variabel, dilakukan secara dinamis.

Tujuan:

Mendapatkan model yang sederhana namun cukup akurat dalam memprediksi nilai respon.
Pemilihan antara pendekatan Konfirmatori atau Eksploratori bergantung pada tujuan penelitian. Jika terdapat teori yang jelas dan mendasar, maka pendekatan konfirmatori lebih sesuai. Namun, jika belum ada teori yang kuat, pendekatan eksploratori bisa membantu menemukan hubungan yang relevan untuk diuji lebih lanjut dalam penelitian selanjutnya.

Simulasi Data

# Membangun Model Regresi Logistik
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
## Warning: package 'lattice' was built under R version 4.3.3
## 
## 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)

## Simulasi Data
set.seed(57)
n <- 200
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
lin_pred <- -0.5 + 1.5 * x1 + 0.4 * x2 - 0.9 * 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
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.5254     0.2708  -1.941   0.0523 .  
## x1            1.4894     0.2315   6.433 1.25e-10 ***
## x2            0.4708     0.3692   1.275   0.2023    
## x3           -0.8609     0.1982  -4.343 1.41e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.28  on 199  degrees of freedom
## Residual deviance: 190.58  on 196  degrees of freedom
## AIC: 198.58
## 
## Number of Fisher Scoring iterations: 5

Metode Stepwise: Forward, Backward, dan Kedua Arah

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

Interpretasi: Model dengan AIC lebih rendah umumnya lebih disukai untuk dipilih menjadi model terbaik.

Evaluasi Model: ROC dan AUC

## 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 = "red")

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

Pseudo R-Squared

## Pseudo R-Squared
PseudoR2(step_both, which = c("CoxSnell", "Nagelkerke", "McFadden"))
##   CoxSnell Nagelkerke   McFadden 
##  0.3431234  0.4582487  0.3042290

Interpretasi: 1. Cox & Snell (0.2386): Nilai ini menunjukkan bahwa model menjelaskan sekitar 23.86% variasi dalam data. Namun, nilai ini memiliki batas atas yang lebih rendah dari 1, sehingga interpretasinya terbatas.

  1. Nagelkerke (0.3294): Merupakan penyesuaian dari Cox & Snell agar skalanya berkisar antara 0 hingga 1. Dengan nilai 32.94%, ini menunjukkan bahwa model menjelaskan proporsi variasi yang cukup baik dan lebih mudah diinterpretasikan dibandingkan Cox & Snell.

  2. McFadden (0.2115): Mengindikasikan bahwa model memiliki kemampuan prediktif yang cukup layak. Nilai antara 0.2 hingga 0.4 umumnya dianggap sebagai indikasi model yang baik dalam konteks regresi logistik.

Tabel Klasifikasi dan Evaluasi

## 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 86 25
##          1 21 68
##                                           
##                Accuracy : 0.77            
##                  95% CI : (0.7054, 0.8264)
##     No Information Rate : 0.535           
##     P-Value [Acc > NIR] : 5.176e-12       
##                                           
##                   Kappa : 0.5364          
##                                           
##  Mcnemar's Test P-Value : 0.6583          
##                                           
##             Sensitivity : 0.7312          
##             Specificity : 0.8037          
##          Pos Pred Value : 0.7640          
##          Neg Pred Value : 0.7748          
##              Prevalence : 0.4650          
##          Detection Rate : 0.3400          
##    Detection Prevalence : 0.4450          
##       Balanced Accuracy : 0.7675          
##                                           
##        'Positive' Class : 1               
## 
conf_matrix$byClass[c("Sensitivity", "Specificity")]
## Sensitivity Specificity 
##   0.7311828   0.8037383

Model terbaik berdasarkan AIC adalah model step both dengan AIC paling minimum. Dengan nilai ROC menunjukan 0.8444 yang terbilang cukup baik. Nilai Pseudo R-Squared dan hasil klasifikasi model juga menunjukkan hasil yang cuku baik. Dimana klasifikasi memiliki akurasi sebesar 77%.

Metode Perbandingan Model dalam Regresi Logistik

Tujuan Dalam sub-bab ini akan disajikan cara membandingkan model regresi logistik menggunakan ukuran Deviance, AIC, dan Likelihood-Ratio serta menjelaskan prinsip Parsimony dalam pemilihan model.

Simulasi Data

# Metode Perbandingan Model dalam Regresi Logistik
library(MASS)
library(broom)
library(DescTools)

## Simulasi Data
set.seed(57)
n <- 300
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
lin_pred <- -1 + 1.5 * x1 - 0.7 * x2 + 0.9 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)

Pembuatan Model

## Pembuatan Model
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)

Perbandingan AIC dan Deviance

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

Dari perbandingan nilai AIC dan Deviance dari 3 model, model terbaik adalah model dengan AIC dan Deviance paling kecil. Hal itu terdapat pada model dengan variabel paling banyak.

Likelihood-Rasio Test

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

Uji Likelihood-Rasio Test menunjukkan p-value 0.001768. P-value < 0.05 maka secara signifikan penambahan variabel x2 meningkatkan kecocokan pada model. x2 layak dipertahankan.

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

Uji Likelihood-Rasio Test menunjukkan p-value 4.849e-11. P-value < 0.05 maka secara signifikan penambahan variabel x2 meningkatkan kecocokan pada model. x3 layak dipertahankan.

Prinsip Parsimony

Dalam pemodelan, model yang kompleks umumnya menghasilkan nilai AIC dan deviance yang lebih kecil. Namun, model yang lebih sederhana lebih mudah untuk dipahami dan dijelaskan. Jika perbedaan AIC antar model tidak signifikan, maka model yang lebih sederhana sebaiknya dipilih. Model kompleks yang terlalu overfit bisa jadi tidak generalis.

Rumus dan Penjabaran

Rumus AIC

\[ \text{AIC} = -2(\log L - k) = -2 \log L + 2k \]

Penjelasan:
AIC (Akaike Information Criterion) digunakan untuk menilai seberapa baik model sesuai dengan data (melalui likelihood) dengan penalti untuk kompleksitas model (jumlah parameter, \(k\)).
Model dengan nilai AIC lebih kecil dianggap lebih baik. Namun, jika perbedaannya tidak signifikan, pilihlah model yang lebih sederhana, karena model yang terlalu kompleks rentan terhadap overfitting.

Rumus Deviance

\[ D = -2 \log(L_{\text{model}}) - \log(L_{\text{model\_saturasi}}) \]

Penjelasan: Deviance adalah ukuran seberapa jauh model saat ini dari model yang sangat sesuai (model saturasi). Semakin kecil nilai deviance, semakin baik model tersebut dalam menjelaskan variasi dalam data.

Rumus Likelihood Ratio Test

\[ G^2 = -2 \log L_0 + 2 \log L_1 \]

Penjelasan: Statistik Likelihood Ratio Test (LRT) digunakan untuk membandingkan dua model yang bersarang (nested). Model yang lebih kompleks diuji apakah memberikan peningkatan signifikan terhadap kecocokan model dibanding model yang lebih sederhana.
Jika nilai \(G^2\) besar dan p-value signifikan, maka model kompleks dianggap lebih baik dari model sederhana.

Evaluasi Tabel dan Klasifikasi dan Akurasi Model

## Evaluasi Tabel Klasifikasi dan Akurasi Model
pred_prob <- predict(model3, type = "response")
pred_class <- factor(ifelse(pred_prob >= 0.5, 1, 0))
conf_matrix <- confusionMatrix(pred_class, data$y, positive = "1")
conf_matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 194  38
##          1  18  50
##                                           
##                Accuracy : 0.8133          
##                  95% CI : (0.7646, 0.8558)
##     No Information Rate : 0.7067          
##     P-Value [Acc > NIR] : 1.589e-05       
##                                           
##                   Kappa : 0.5177          
##                                           
##  Mcnemar's Test P-Value : 0.01112         
##                                           
##             Sensitivity : 0.5682          
##             Specificity : 0.9151          
##          Pos Pred Value : 0.7353          
##          Neg Pred Value : 0.8362          
##              Prevalence : 0.2933          
##          Detection Rate : 0.1667          
##    Detection Prevalence : 0.2267          
##       Balanced Accuracy : 0.7416          
##                                           
##        'Positive' Class : 1               
## 

Sensitivitas dan Spesifisitas

  • Sensitivitas: Kemampuan model mendeteksi kelas positif secara benar (True Positive Rate) \[ \text{Sensitivity}=\frac{TP}{TP+FN} \] \[ \text{Sensitivity}=\frac{50}{50+38}=0.568 \]
  • Spesifisitas: Kemampuan model mendeteksi kelas negatif secara benar (True Negative Rate) \[ \text{Spesifisitas}=\frac{TN}{TN+FP} \] \[ \text{Spesifisitas}=\frac{194}{194+18}=0.915 \]
conf_matrix$byClass[c("Sensitivity", "Specificity")]
## Sensitivity Specificity 
##   0.5681818   0.9150943

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

Detail ROC, Penjelasan Kurva ROC (Receiver Operating Characteristic)

  1. Pengertian Kurva ROC

Kurva ROC (Receiver Operating Characteristic) merupakan alat visual yang sering digunakan untuk mengevaluasi performa dari model klasifikasi biner. Kurva ini memperlihatkan hubungan antara tingkat True Positive Rate (Sensitivity) dan False Positive Rate (1 - Specificity) pada berbagai nilai ambang klasifikasi.

  • Sumbu Y (Sensitivity): TP / (TP + FN)
  • Sumbu X (1 - Specificity): FP / (FP + TN)
  • Garis diagonal (dari kiri bawah ke kanan atas) menggambarkan model tanpa kemampuan klasifikasi (seperti tebakan acak).
  • Model yang ideal akan memiliki kurva mendekati sudut kiri atas, menandakan performa klasifikasi yang baik.
  1. Nilai Ambang dan Pergeseran Kurva

Perubahan nilai ambang klasifikasi (cut-off) dapat mempengaruhi banyaknya data yang diklasifikasikan sebagai positif:

  • Ketika cut-off diturunkan:
    • Sensitivitas meningkat
    • Spesifisitas menurun
  • Ketika cut-off dinaikkan:
    • Sensitivitas menurun
    • Spesifisitas meningkat
  1. Karakteristik Kurva ROC Ideal

Kurva ROC yang ideal memiliki ciri-ciri berikut:

  • Naik secara vertikal tajam hingga sensitivitas mencapai 1
  • Bergerak horizontal ke arah kanan hingga 1 - spesifisitas = 1
  • Area di bawah kurva (AUC) mendekati 1
  1. Interpretasi Nilai AUC (Area Under the Curve)

Interpretasi dari nilai AUC adalah sebagai berikut:

  • AUC = 0.5: Model tidak lebih baik dari tebakan acak
  • AUC > 0.7: Model tergolong cukup baik
  • AUC > 0.9: Model sangat baik

Nilai AUC juga bisa diinterpretasikan sebagai probabilitas bahwa model dapat memberikan skor yang lebih tinggi kepada kasus positif dibandingkan kasus negatif.

  1. Kegunaan Kurva ROC

Kurva ROC digunakan untuk:

  • Membandingkan performa dari beberapa model klasifikasi
  • Menentukan nilai ambang klasifikasi (threshold) terbaik sesuai kebutuhan spesifik (misalnya, ketika false negative lebih berbahaya dibanding false positive)
  1. Visualisasi di R

Di dalam R, kurva ROC dapat dibuat menggunakan package pROC.

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

auc(roc_obj)
## Area under the curve: 0.8389
  1. Simulasi Pemilihan Treshold Optimal Untuk dapat memilih treshold optimal/terbaik, bisa dilakukan dengan mengevaluasi sensitivitas dan spesifisitas pada berbagai cut-off
thresholds <- seq(0.1, 0.9, by = 0.05)
results <- data.frame(Threshold = thresholds)

results$Sensitivity <- sapply(thresholds, function(t) {
pred_class <- ifelse(pred >= t, 1, 0)
cm <- table(Pred = pred_class, Obs = data$y)
TP <- cm["1", "1"]
FN <- cm["0", "1"]
TP / (TP + FN)
})
results$Specificity <- sapply(thresholds, function(t) {
pred_class <- ifelse(pred >= t, 1, 0)
cm <- table(Pred = pred_class, Obs = data$y)
TN <- cm["0", "0"]
FP <- cm["1", "0"]
TN / (TN + FP)
})
print(results)
##    Threshold Sensitivity Specificity
## 1       0.10   0.9866667       0.336
## 2       0.15   0.9200000       0.448
## 3       0.20   0.8800000       0.504
## 4       0.25   0.8133333       0.600
## 5       0.30   0.8133333       0.672
## 6       0.35   0.8000000       0.704
## 7       0.40   0.7333333       0.736
## 8       0.45   0.7200000       0.824
## 9       0.50   0.6800000       0.864
## 10      0.55   0.6133333       0.888
## 11      0.60   0.5733333       0.904
## 12      0.65   0.5066667       0.928
## 13      0.70   0.4000000       0.984
## 14      0.75   0.3200000       0.992
## 15      0.80   0.2666667       0.992
## 16      0.85   0.1866667       1.000
## 17      0.90   0.0800000       1.000

Cut-off optimal dapat dipilih berdasarkan: - Nilai maksimum dari Sensitivity + Spesifisitas - Mempertimbangkan trade-off sesuai tujuan aplikasi

Penggunaan ROC cocok digunakan saat proporsi kelas seimbang. Namun, jika tidak precision-recall curve bisa lebih informatif.

Precision-Recall Curve (PR Curve)

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

  1. Pengertian
  • Presisi (Precision): Merupakan rasio prediksi positif yang benar, yaitu berapa banyak dari seluruh prediksi positif yang benar-benar relevan. \[ \text{Precision} = \frac{TP}{TP + FP} \]

  • Recall (Sensitivitas): Merupakan rasio kasus positif yang berhasil diidentifikasi dengan benar sebagai positif. \[ \text{Recall} = \frac{TP}{TP + FN} \]

  1. Interpretasi
  • Kurva Precision-Recall (PR) menunjukkan perubahan nilai presisi ketika recall meningkat.
  • Tujuan ideal adalah memiliki presisi dan recall yang tinggi, meskipun dalam praktik sering terjadi kompromi (trade-off) antara keduanya.
  • Model yang baik ditunjukkan dengan kurva PR yang cenderung mendekati sudut kanan atas grafik.
  1. Luas Area di Bawah Kurva PR (AUPRC)
  • Area di bawah kurva PR (AUPRC) mendekati nilai 1 menunjukkan performa model yang sangat baik.
  • Nilai baseline AUPRC umumnya setara dengan proporsi kelas positif dalam dataset (prevalensi).
  1. Perbandingan Kurva PR dan Kurva ROC
Aspek Kurva ROC Kurva Precision-Recall
Fokus Semua kelas Kelas positif saja
Kuat di Data seimbang Data tidak seimbang
Sumbu Y Sensitivitas (Recall) Precision
Sumbu X 1 - Spesifisitas Recall
  1. Visualisasi PR Curve di R

Kurva Precision-Recall dapat divisualisasikan menggunakan berbagai package di R, seperti PRROC atau precrec, untuk menganalisis performa model dalam konteks ketidakseimbangan data.

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(57)
x1 <- rnorm(200)
x2 <- rbinom(200, 1, 0.5)
x3 <- rnorm(200)
lin_pred <- -1 + 1.3 * x1 - 0.5 * x2 + 0.7 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(200, 1, p)
data <- data.frame(y = y, x1, x2, x3)
model <- glm(y ~ x1 + x2 + x3, data = data, family = binomial)
prob <- predict(model, type = "response")
pr <- pr.curve(scores.class0 = prob[data$y == 1],
scores.class1 = prob[data$y == 0],
curve = TRUE)
plot(pr)

Pr Curve sangat informatif untuk aplikasi seperti deteksi penipuan atau diagnosis penyakit langka. PR Curve cocok digunakan saat: - Kelas positiv jauh lebih jarang dibanding kelas negatif - Tujuan aplikasi lebih mementingkan presisi terhadap kelas minoritas

Pseudo R-Squared pada Regresi Logistik

Tujuan Pada sub-bab ini akan dijelaskan dan dihitung pseudo R-squared dalam regresi logistik: - \(R_{CoxandSnell}^2\) - \(R_{McFadden}^2\)

Simulasi Data

set.seed(57)
n <- 300
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)
lin_pred <- -1 + 1.4 * x1 - 0.4 * 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_{CoxandSnell}^2=1-(\frac{L_0}{L_M})^{2/n} \] \[ R_{McFadden}^2=1-\frac{logL_M}{logL_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 R Menggunakan pscl

library(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 
## -124.6653690 -187.2993920  125.2680460    0.3344059    0.3413481    0.4786741

Menggunakan DescTools

if (!require(DescTools)) install.packages("DescTools"); library(DescTools)
PseudoR2(model, which = "all")
##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.3344059       0.3130497       0.3413481       0.4786741       0.2945626 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.5304650       0.3769812       0.5694989       0.3775846     257.3307379 
##             BIC          logLik         logLik0              G2 
##     272.1458678    -124.6653690    -187.2993920     125.2680460

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

Apa itu Distribusi Multinomial

  1. Pengantar

Distribusi multinomial merupakan generalisasi dari distribusi binomial untuk kasus dengan lebih dari dua kategori. Jika suatu percobaan menghasilkan k kategori, dan kita ingin mengetahui probabilitas dari sejumlah kejadian dalam setiap kategori, maka distribusi multinomial bisa digunakan.

Jika variabel acak \(( X_1, X_2, ..., X_k)\) menyatakan jumlah kejadian pada masing-masing dari \(( k )\) kategori, maka distribusi peluangnya diberikan oleh:

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

dengan syarat:

  1. Studi Kasus

Sebuah survei dilakukan terhadap 10 responden yang diminta memilih satu dari tiga buah favorit: Apel (A), Jeruk (J), dan Pisang (P).

Hasil Survei: - Apel: 3 orang
- Jeruk: 5 orang
- Pisang: 2 orang

Probabilitas Teoritis Preferensi: - \(( p_A = 0.3 )\)
- \(( p_J = 0.5 )\)
- \(( p_P = 0.2 )\)

Pertanyaan: Berapa peluang bahwa 10 orang akan terdiri dari 3 pemilih Apel, 5 pemilih Jeruk, dan 2 pemilih Pisang?

Perhitungan Manual di R

n <- 10
x <- c(3, 5, 2)
p <- c(0.3, 0.5, 0.2)
# 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.08505

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

Regresi Logistik Multinomial

Regresi logistik multinomial digunakan ketika kita ingin menganalisis hubungan antara variabel respon kategorik yang memiliki lebih dari dua kategori dengan satu atau lebih variabel prediktor.

Misalkan variabel respon \(( Y )\) memiliki \(( K )\) kategori, dan kita memilih salah satu kategori (biasanya kategori terakhir) sebagai acuan (baseline), maka model logit untuk kategori \(( j )\) terhadap kategori baseline \(( K )\) adalah:

\[ \log\left(\frac{P(Y = j)}{P(Y = K)}\right) = \beta_{0j} + \beta_{1j}x_1 + \cdots + \beta_{pj}x_p \]

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

Model Logit Baseline-Kategori

Model baseline-category logit membandingkan log odds (logit) dari setiap kategori terhadap satu kategori acuan (baseline). Bentuk umumnya adalah:

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

di mana:

  • \(( \pi_j )\) adalah probabilitas kategori ke-j
  • \(( \pi_c )\) adalah probabilitas kategori referensi (baseline)

Dengan demikian, akan terbentuk sebanyak \(( (c - 1) )\) fungsi logit.

Catatan: Dalam R, kategori referensi biasanya otomatis ditetapkan sebagai kategori terakhir, namun bisa diubah menggunakan fungsi relevel().

Model Regresi untuk Satu Prediktor

Jika hanya terdapat satu prediktor \(( x )\), maka bentuk persamaan regresi logit menjadi:

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

Contoh Kasus: Respon dengan 3 Kategori

Misalnya variabel \(( Y )\) memiliki tiga kategori: \(( Y \in \{1, 2, 3\} )\), dan kategori ke-3 dijadikan acuan. Maka modelnya:

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

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

Relasi Antar Kategori Untuk melihat hubungan antara kategori 1 dan 2, kita bisa menghitung logit relatifnya:

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

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

Dengan ini, kita dapat membandingkan logit antara dua kategori non-baseline.

Implementasi di R Model ini dapat diimplementasikan dengan fungsi multinom() dari package nnet.

Estimasi Parameter

Parameter dalam model ini biasanya diestimasi menggunakan metode Maximum Likelihood, dengan pendekatan numerik seperti algoritma Newton-Raphson.

Fungsi log-likelihood untuk data dengan \(( n )\) observasi dan \(( k )\) kategori:

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

dengan:

  • \(( \pi_{ij} = P(Y_i = j \mid x_i) )\)
  • \(( y_{ij} = 1 )\) jika observasi ke-\(( i )\) berada pada kategori ke-\(( j )\), dan 0 selainnya.

Contoh Kasus

Sebuah restoran cepat saji ingin mengetahui faktor-faktor yang memengaruhi preferensi pelanggan terhadap jenis minuman yang dipesan, yaitu: Teh, Kopi, atau Jus.

Restoran melakukan survei terhadap 200 pelanggan dan mengumpulkan data berikut:

  • Drink: Jenis minuman yang dipesan (Teh, Kopi, Jus)
  • Age: Usia pelanggan
  • TimeOfDay: Waktu kunjungan (Pagi, Siang, Malam)
  • VisitFrequency: Frekuensi kunjungan dalam sebulan

Tujuannya adalah: Mengetahui bagaimana usia, waktu kunjungan, dan frekuensi kunjungan memengaruhi pilihan jenis minuman pelanggan.

Simulasi Data

set.seed(57)
n <- 200
TimeofDay <- sample(c("Pagi", "Siang", "Malam"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 30, sd = 8))
VisitFrequency <- round(pmax(rnorm(n, mean = 10, sd = 3), 0))
# Simulasikan Time of Day berdasarkan probabilitas berbeda per Jenis Minuman
Drink <- sapply(TimeofDay, function(time) {
if (time == "Pagi") {
sample(c("Teh", "Kopi", "Jus"), size = 1, prob = c(0.3, 0.4, 0.3))
} else if (time == "Siang") {
sample(c("Teh", "Kopi", "Jus"), size = 1, prob = c(0.2, 0.4, 0.4))
} else {
sample(c("Teh", "Kopi", "Jus"), size = 1, prob = c(0.4, 0.4, 0.2))
}
})
df <- data.frame(Drink = factor(Drink), Age, TimeofDay = factor(TimeofDay), VisitFrequency)
df$Drink <- relevel(df$Drink, ref = "Teh") # baseline
head(df)

Estimasi Model

library(nnet)
## Warning: package 'nnet' was built under R version 4.3.3
model_mnlogit <- multinom(Drink ~ Age + TimeofDay + VisitFrequency, data = df)
## # weights:  18 (10 variable)
## initial  value 219.722458 
## iter  10 value 209.542916
## final  value 209.300869 
## converged
summary(model_mnlogit)
## Call:
## multinom(formula = Drink ~ Age + TimeofDay + VisitFrequency, 
##     data = df)
## 
## Coefficients:
##      (Intercept)          Age TimeofDayPagi TimeofDaySiang VisitFrequency
## Jus  -1.26942303 0.0074754230     1.0855274      1.4624395    0.027840437
## Kopi -0.01732173 0.0004726801     0.7983403      0.8789951   -0.003789478
## 
## Std. Errors:
##      (Intercept)        Age TimeofDayPagi TimeofDaySiang VisitFrequency
## Jus    0.9710145 0.02200367     0.4878910      0.4908274     0.06270705
## Kopi   0.8347962 0.01969198     0.4224498      0.4376284     0.05616052
## 
## Residual Deviance: 418.6017 
## AIC: 438.6017

Nilai P-Value dan Interpretasi

z <- summary(model_mnlogit)$coefficients / summary(model_mnlogit)$standard.errors
pval <- 2 * (1 - pnorm(abs(z)))
round(pval, 4)
##      (Intercept)    Age TimeofDayPagi TimeofDaySiang VisitFrequency
## Jus       0.1911 0.7341        0.0261         0.0029         0.6571
## Kopi      0.9834 0.9808        0.0588         0.0446         0.9462

Interpretasi: Koefisien untuk kategori “Kopi” dan “Jus” dibandingkan dengan baseline “Teh”. Jika nilai p-value < 0.05 menunjukkan bahwa variabel tersebut signifikan memengaruhi preferensi minuman yang dipesan di restauran tersebut.

Prediksi dan Validasi

df$Predicted <- predict(model_mnlogit)
table(Predicted = df$Predicted, Actual = df$Drink)
##          Actual
## Predicted Teh Jus Kopi
##      Teh   26  12   25
##      Jus    1   2    4
##      Kopi  27  45   58

Kesimpulan

Model regresi logistik multinomial terbukti efektif dalam: - Mengkaji keterkaitan antara karakteristik pelanggan dan preferensi pemilihan jenis minuman

  • Mengidentifikasi faktor-faktor yang secara signifikan memengaruhi keputusan pemilihan

  • Memberikan kemampuan untuk memprediksi jenis minuman yang akan dipesan oleh pelanggan berdasarkan atribut pribadi mereka

Regresi Logistik Ordinal

Regresi logistik ordinal digunakan saat variabel respons memiliki sifat urutan, contohnya seperti tingkat kepuasan pelanggan: Tidak Puas, Cukup Puas, dan Sangat Puas.

Model ini berbeda dengan:

Konsep Model Cumulative Logit

Model yang digunakan adalah Cumulative Logit Model dengan asumsi proportional odds. Model ini memodelkan probabilitas kumulatif hingga kategori tertentu:

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

Penjelasan: - \(\alpha_j\): intercept khusus untuk batas kategori ke-\(j\) - \(\beta\): koefisien regresi (diasumsikan sama untuk semua batas kategori)

Jika ada \(k\) kategori, maka dibentuk \((k - 1)\) model logit kumulatif.

Interpretasi Koefisien

Koefisien \(\beta\) menunjukkan arah dan besar pengaruh dari variabel prediktor \(x\) terhadap peluang berada pada kategori lebih rendah atau sama.

  • Jika \(\beta > 0\): kenaikan \(x\) meningkatkan kemungkinan individu berada pada kategori lebih rendah
  • Jika \(\beta < 0\): kenaikan \(x\) meningkatkan kemungkinan individu berada pada kategori lebih tinggi

Odds ratio (rasio peluang) dihitung sebagai:

\[ \text{OR} = e^{\beta} \] ## Contoh Kasus Sebuah kampus swasta melakukan survei terhadap mahasiswanya untuk mengevaluasi tingkat kepuasan terhadap fasilitas kampus. Tujuan dari survei ini adalah untuk mengetahui apakah jumlah fasilitas yang tersedia di kampus (seperti ruang belajar, Wi-Fi, kantin, laboratorium, area olahraga, dll.) berkaitan dengan kepuasan mahasiswa terhadap layanan dan kenyamanan lingkungan kampus.

Dalam penelitian ini, kampus dikelompokkan berdasarkan jumlah fasilitas yang tersedia menjadi tiga kategori: - Sedikit (1-3 fasilitas) - Sedang (4-6 fasilitas) - Banyak (7 fasilitas atau lebih)

Sementara itu, tingkat kepuasan mahasiswa dikategorikan menjadi tiga tingkat: - Tidak puas - Cukup - Puas

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

Estimasi Model Ordinal

model_ord <- polr(satisfaction ~ quality, data = df, Hess = TRUE)
summary(model_ord)
## Call:
## polr(formula = satisfaction ~ quality, data = df, Hess = TRUE)
## 
## Coefficients:
##          Value Std. Error t value
## quality 0.5442    0.06487   8.389
## 
## Intercepts:
##                  Value  Std. Error t value
## Tidak Puas|Cukup 1.9201 0.3488     5.5053 
## Cukup|Puas       4.6714 0.4916     9.5033 
## 
## Residual Deviance: 336.6941 
## AIC: 342.6941

Nilai P-Value

(ctable <- coef(summary(model_ord)))
##                      Value Std. Error  t value
## quality          0.5442158 0.06486982 8.389352
## Tidak Puas|Cukup 1.9200523 0.34876696 5.505259
## Cukup|Puas       4.6713599 0.49155106 9.503305
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
## quality          0.5442158 0.06486982 8.389352       0
## Tidak Puas|Cukup 1.9200523 0.34876696 5.505259       0
## Cukup|Puas       4.6713599 0.49155106 9.503305       0

Prediksi Probabilitas

newdata <- data.frame(quality = 5:9)
predict(model_ord, newdata = newdata, type = "probs")
##   Tidak Puas     Cukup      Puas
## 1 0.30980594 0.5656713 0.1245227
## 2 0.20664925 0.5964953 0.1968554
## 3 0.13130624 0.5717406 0.2969531
## 4 0.08064066 0.4981067 0.4212526
## 5 0.04843479 0.3951607 0.5564045

Goodness-of-Fit dan Propotional 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 model cumulative logit, terdapat beberapa pendekatan lain untuk menangani data ordinal:

  • Adjacent-category logit: membandingkan kategori satu dengan yang berdekatan
  • Continuation-ratio logit: melihat kemungkinan pindah ke kategori berikutnya secara berurutan

Model alternatif ini bisa digunakan apabila asumsi proportional odds tidak terpenuhi pada cumulative logit.

Kesimpulan

  • Regresi ordinal cocok digunakan ketika variabel respons memiliki urutan.
  • Model cumulative logit memodelkan probabilitas kumulatif dalam bentuk log-odds.
  • Di R, model ini bisa diimplementasikan dengan fungsi polr() dari paket MASS.
  • Jika ingin melakukan validasi model, bisa digunakan uji devians atau likelihood ratio test.

Asumsi Paralelisme dalam Regresi Logistik Ordinal

Regresi logistik ordinal yang paling umum digunakan adalah Cumulative Logit Model dengan asumsi proportional odds.

Asumsi ini juga dikenal sebagai asumsi paralelisme (parallel lines assumption), yang menyatakan bahwa: Koefisien regresi \(( \beta )\) tetap sama untuk setiap batas kumulatif dari kategori respons.

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

Ciri-cirinya: - Intercept \(( \alpha_j )\) berbeda untuk setiap batas kategori - Koefisien \(( \beta )\) tetap, menunjukkan efek variabel prediktor sama untuk semua kategori

Visualisasi asumsi paralelisme: - Jika asumsi ini terpenuhi, garis logit kumulatif terhadap \(( x )\) akan sejajar (slope sama, hanya beda posisi)

Jika Asumsi Tidak Terpenuhi: - Koefisien regresi bisa bervariasi di setiap batas kategori - Gunakan model yang lebih fleksibel, seperti: - Partial Proportional Odds Model - Generalized Ordinal Logistic Regression

Pengujian Asumsi Paralelisme

Untuk menguji apakah asumsi paralelisme terpenuhi, dapat digunakan:

  • Likelihood Ratio Test: membandingkan model proportional vs non-proportional
  • Brant Test (paket brant di R)

Contoh dengan package brant di R:

library(MASS)
library(brant)
## Warning: package 'brant' was built under R version 4.3.3
brant(model_ord)
## -------------------------------------------- 
## Test for X2  df  probability 
## -------------------------------------------- 
## Omnibus      0.46    1   0.5
## quality      0.46    1   0.5
## -------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds

Jika p-value < 0.05 maka asumsi nya 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 berperan penting dalam statistik terapan karena banyak fenomena di dunia nyata menghasilkan data dalam bentuk kategori, seperti jenis kelamin, status pekerjaan, dan tingkat pendidikan.

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

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(40, 15, 35, 60), nrow=2, 
       dimnames = list(Diet = c("Rendah Karbo", "Rendah Lemak"),
                       Hasil = c("Ya", "Tidak")))
table_data
##               Hasil
## Diet           Ya Tidak
##   Rendah Karbo 40    35
##   Rendah Lemak 15    60

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

library(MASS)
loglm(~ Diet * Hasil, data = table_data)
## Call:
## loglm(formula = ~Diet * Hasil, 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 \]

data_glm <- data.frame(
  Hasil = c(1, 0, 1, 0),
  Diet = factor(c("Rendah Karbo", "Rendah Karbo", "Rendah Lemak", "Rendah Lemak")),
  Frek = c(40, 15, 35, 60)
)
model_logit <- glm(Hasil ~ Diet, weights = Frek, family = binomial, data = data_glm)
summary(model_logit)
## 
## Call:
## glm(formula = Hasil ~ Diet, family = binomial, data = data_glm, 
##     weights = Frek)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        0.9808     0.3028   3.240   0.0012 ** 
## DietRendah Lemak  -1.5198     0.3700  -4.108    4e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 207.94  on 3  degrees of freedom
## Residual deviance: 189.50  on 2  degrees of freedom
## AIC: 193.5
## 
## 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(40, 15, 35, 60), nrow=2,
       dimnames = list(Diet = c("Rendah Karbo", "Rendah Lemak"),
                       Hasil = c("Ya", "Tidak")))
##               Hasil
## Diet           Ya Tidak
##   Rendah Karbo 40    35
##   Rendah Lemak 15    60

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(45, 55, 65, 35), nrow=2, byrow=TRUE)
dimnames(data) <- list(Diet = c("Rendah Karbo", "Rendah Lemak"), Hasil = c("Ya", "Tidak"))
ftable(data)
##              Hasil Ya Tidak
## Diet                       
## Rendah Karbo       45    55
## Rendah Lemak       65    35

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

model_saturated <- loglm(~ Diet * Hasil, data = data)
summary(model_saturated)
## Formula:
## ~Diet * Hasil
## attr(,"variables")
## list(Diet, Hasil)
## attr(,"factors")
##       Diet Hasil Diet:Hasil
## Diet     1     0          1
## Hasil    0     1          1
## attr(,"term.labels")
## [1] "Diet"       "Hasil"      "Diet:Hasil"
## 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(~ Diet + Hasil, data = data)
summary(model_indep)
## Formula:
## ~Diet + Hasil
## attr(,"variables")
## list(Diet, Hasil)
## attr(,"factors")
##       Diet Hasil
## Diet     1     0
## Hasil    0     1
## attr(,"term.labels")
## [1] "Diet"  "Hasil"
## 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 8.138435  1 0.004333667
## Pearson          8.080808  1 0.004473649

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 antara 2 kategori
  • OR > 1: Asosiasi positif antara 2 kategori
  • OR < 1: Asosiasi negatif antara 2 kategori

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

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:
##  ~Diet + Hasil 
## Model 2:
##  ~Diet * Hasil 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   8.138435  1                                    
## Model 2   0.000000  0   8.138435         1        0.00433
## Saturated 0.000000  0   0.000000         0        1.00000

Model Log-Linear pada tabel kontingensi

Model log-linear merupakan metode yang digunakan untuk menganalisis hubungan antara dua atau lebih variabel kategorik yang disusun dalam bentuk tabel kontingensi. Model ini bekerja dengan mengasumsikan bahwa logaritma dari nilai harapan frekuensi setiap sel (μᵢⱼ) dapat dinyatakan sebagai penjumlahan dari efek masing-masing variabel serta (jika diperlukan) efek interaksi antar variabel. Untuk tabel 2x2, bentuk modelnya adalah: \[ \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 (jumlah kejadian) dalam tabel kontingensi serta menguji hubungan antar variabel kategorik, tanpa menganggap adanya variabel respon dan prediktor.

  • Model regresi logistik bertujuan untuk memodelkan probabilitas terjadinya suatu outcome (biner) berdasarkan satu atau lebih prediktor, baik yang bersifat kategorik maupun kontinu.

Estimasi Parameter Log Linear 2 Arah

Sistem Persamaan Model Log-Linear \[\begin{align*} \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} \\ \end{align*}\] Constraint Sum-to-Zero \[\begin{align*} \lambda_1^A + \lambda_2^A = 0 \\ \\ \lambda_1^B + \lambda_2^B = 0 \\ \\ \lambda_{11}^{AB} + \lambda_{12}^{AB} + \lambda_{21}^{AB} + \lambda_{22}^{AB} = 0 \end{align*}\] Rumus Estimasi Parameter dengan Sum-to-Zero Constraint

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

Contoh Kasus:

Seorang ilmuwan gizi melakukan penelitian untuk mengetahui apakah diet berpengaruh signifikan terhadap keberhasilan seseorang dalam menurunkan berat badan. Penelitian ini dilakukan terhadap 90 partisipan yang memiliki tujuan menurunkan berat badan dalam waktu 3 bulan.

Partisipan dibagi ke dalam dua kelompok berdasarkan kebiasaan mereka menjalankan program diet: - Melakukan Diet - Tidak Melakukan Diet

Setelah 3 bulan, dicatat apakah masing-masing partisipan berhasil menurunkan berat badan secara signifikan (misalnya, minimal 5% dari berat badan awal) atau tidak berhasil.

Data kemudian dikumpulkan dan disajikan dalam tabel kontingensi 2x2 berikut:

Diet Berhasil Tidak
Ya 40 15
Tidak 10 25

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 \] ## Estimasi Parameter Model (Manual, Sum-to-zero) Misalkan: - A1 = Diet (Ya), A2 = Tidak - B1 = Berhasil, B2 = Tidak Observasi: - \(n_{11}\) = 40, \(n_{12}\) = 15 - \(n_{21}\) = 10, \(n_{22}\) = 25

\[\begin{align*} \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} \\ \end{align*}\]

Constraint sum-to-zero:

\[\begin{align*} \lambda_1^A + \lambda_2^A = 0 \\ \\ \lambda_1^B + \lambda_2^B = 0 \\ \\ \lambda_{11}^{AB} + \lambda_{12}^{AB} + \lambda_{21}^{AB} + \lambda_{22}^{AB} = 0 \end{align*}\]

Langkah-langkah: 1. Hitung rata-rata log frekuensi sel: \[\begin{align*} \lambda &= \frac{1}{4} \sum_{i=1}^{2} \sum_{j=1}^{2} \log(n_{ij}) \\ &= \frac{1}{4} [\log(40) + \log(15) + \log(10) + \log(25)] \\ &= 2.979 \end{align*}\]

  1. Efek Utama A (Merokok): \[\begin{align*} \lambda_1^A &= \frac{1}{2} [(\log(40) + \log(15)) - (\log(10) + \log(25))] \\ \\ &= \frac{1}{2} [(3.6889 + 2.7081) - (2.3026 + 3.2189)] \\ \\ &= \frac{1}{2} (6.3970 - 5.5215) \\ \\ &= \frac{1}{2} (0.8755) \\ \\ &= 0.4378 \\ \\ \lambda_2^A &= -0.4378 \\ \end{align*}\]

  2. Efek Utama B (Status): \[\begin{align*} \\ \lambda_1^B &= \frac{1}{2} [(\log(40) + \log(10)) - (\log(15) + \log(25))] \\ \\ &= \frac{1}{2} [(3.6889 + 2.3026) - (2.7081 + 3.2189)] \\ \\ &= \frac{1}{2} (5.9915 - 5.9270) \\ \\ &= \frac{1}{2} (0.0645) \\ \\ &= 0.0323 \\ \\ \lambda_2^B &= -0.0323 \\ \end{align*}\]

  3. Efek Interaksi: \[\begin{align*} \lambda_{11}^{AB} &= \frac{1}{4} [(\log(40) - \log(15)) - \log(10) + \log(25)] \\ \\ &= \frac{1}{4} [(3.6889 - 2.7081 - 2.3026 + 3.2189)] \\ \\ &= \frac{1}{4} (-1.3218 + 3.2189) \\ \\ &= \frac{1}{4} (1.8971) \\ \\ &= 0.4743 \\ \\ \lambda_{12}^{AB} &= -0.4743 \\ \\ \lambda_{21}^{AB} &= -0.4743 \\ \\ \lambda_{22}^{AB} &= 0.4743 \\ \end{align*}\]

Ringkasan Parameter

  • \(\lambda = 2.979\)
  • \(\lambda_1^A = 0.4378, \lambda_2^A = -0.4378\)
  • \(\lambda_1^B = 0.0323, \lambda_2^B = -0.0323\)
  • \(\lambda_{11}^{AB} = 0.4743, \lambda_{12}^{AB} = -0.4743, \lambda_{21}^{AB} = -0.4743, \lambda_{22}^{AB} = 0.4743\)

Hitung Odds Ratio dan Interval Kepercayaan

\[ OR = \frac{n_{11} n_{22}}{n_{12} n_{21}} = \frac{40 \times 25}{15 \times 10} = \frac{1000}{150} = 6.67 \] Log odds ratio: \[ \log(OR) = \log(6.67) = 1.8971 \] Standard error (SE):

\[ SE = \sqrt{ \frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} } \\ \\ = \sqrt{ \frac{1}{40} + \frac{1}{15} + \frac{1}{10} + \frac{1}{25} }\\ \\ = \sqrt{0.025 + 0.067 + 0.1 + 0.04}\\ \\ = \sqrt{0.232}\\ \\ = 0.4817 \\ \]

95% Confidence Interval for log(OR):

\[\begin{align*} \\ \log(OR) \pm 1.96 \times SE &= 1.8971 \pm 1.96 \times 0.4817 \\ \\ &= 1.8971 \pm 0.9441 \\ \\ &= (1.8971 - 0.9441,\ 1.8971 + 0.9441) \\ \\ &= (0.9530,\ 2.8412) \\ \end{align*}\]

Back-transform to get CI for OR:

\[\begin{align*} \\ Lower &= \exp(0.9535) = 2.5935\\ \\ Upper &= \exp(2.8412) = 17.136 \\ \end{align*}\]

Jadi, OR=6.67 dengan taraf signifikansi 95% memiliki interval kepercayaan dalam rentang 2.5935 - 17.136

Fitting Model Log-Linear dengan R

# Data 2x2
tabel <- matrix(c(40, 15, 10, 25), nrow = 2, byrow = TRUE)
colnames(tabel) <- c("Berhasil", "Tidak")
rownames(tabel) <- c("Ya", "Tidak")
tabel
##       Berhasil Tidak
## Ya          40    15
## Tidak       10    25
#Ubah menjadi bentuk tabel untuk glm
data <- as.data.frame(as.table(tabel))
colnames(data) <- c("Diet", "Hasil", "Freq")
data
# Model tanpa interaksi
fit_no_inter <- glm(Freq ~ Diet + Hasil, family = poisson, data = data)
summary(fit_no_inter)
## 
## Call:
## glm(formula = Freq ~ Diet + Hasil, family = poisson, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.4195     0.1645  20.784   <2e-16 ***
## DietTidak    -0.4520     0.2162  -2.090   0.0366 *  
## HasilTidak   -0.2231     0.2121  -1.052   0.2928    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 22.915  on 3  degrees of freedom
## Residual deviance: 17.319  on 1  degrees of freedom
## AIC: 42.628
## 
## Number of Fisher Scoring iterations: 4
# Model dengan interaksi
fit_inter <- glm(Freq ~ Diet * Hasil, family = poisson, data = data)
summary(fit_inter)
## 
## Call:
## glm(formula = Freq ~ Diet * Hasil, family = poisson, data = data)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            3.6889     0.1581  23.331  < 2e-16 ***
## DietTidak             -1.3863     0.3536  -3.921 8.82e-05 ***
## HasilTidak            -0.9808     0.3028  -3.240   0.0012 ** 
## DietTidak:HasilTidak   1.8971     0.4813   3.942 8.10e-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.2915e+01  on 3  degrees of freedom
## Residual deviance: 4.4409e-16  on 0  degrees of freedom
## AIC: 27.309
## 
## 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.8971\) 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):

Rendah Sedang Tinggi
Pekerja 10 25 15
Pengangguran 20 18 5

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: Status Pekerjaan (i=1:Pekerja, i=2:Pengangguran) - B: Tingkat Kepuasan Hidup (j=1:Rendah, j=2:Sedang, j=3:Tinggi) - 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:

\[\begin{align*} \\ \log(\mu_{ij}) &= \lambda \\ \\ &+ \lambda^A_1(Pekerja), \lambda_2^A(Pengangguran)\\ &+ \lambda^B_1(Rendah), \lambda^B_2(Sedang), \lambda^B_3(Tiggi)\\ \\ &+ \lambda^{AB}_{ij}(interaksi jika ada) \\ \end{align*}\]

Langkah-langkah perhitungan:

  1. Hitung rata-rata log frekuensi sel: \[\begin{align*} \\ \lambda &= \frac{1}{6} \sum_{i=1}^{2} \sum_{j=1}^{3} \log(n_{ij}) \\ \\ &= \frac{1}{6} [\log(10) + \log(25) + \log(15) + \log(20) + \log(18) + \log(5)] \\ \\ &= \frac{1}{6} [2.3026 + 3.2189 + 2.7081 + 2.9957 + 2.8904 + 0.6931] \\ \\ &= \frac{1}{6} (14.8088) \\ \\ &= 2.4681 \\ \end{align*}\]

  2. Efek utama A (Jenis Kelamin): \[\begin{align*} \\ \lambda_1^A &= \frac{1}{3} \left[ (\log 10 + \log 25 + \log 15) - (\log 20 + \log 18 + \log 5) \right] \\ \\ &= \frac{1}{3} [ (2.3026 + 3.2189 + 2.7081) - (2.9957 + 2.8904 + 0.6931) ] \\ \\ &= \frac{1}{3} (8.2296 - 6.5792) \\ \\ &= \frac{1}{3} (1.6504) \\ \\ &= 0.5501 \\ \\ \lambda_2^A &= -\lambda_1^A = -0.5501 \\ \end{align*}\]

  3. Efek utama B (Kategori BMI): \[\begin{align*} \\ \lambda_1^B &= \frac{1}{2} \left[ (\log 10 + \log 20) - (\log 25 + \log 18 + \log 15 + \log 5)/2 \right] \\ \\ &= \frac{1}{2} (2.3026 + 2.9957 - 5.2134) \\ \\ &= \frac{1}{2} (0.0849) \\ \\ &= 0.0425 \\ \\ \lambda_2^B &= \frac{1}{2} (3.2189 + 2.8904 - 4.8079) = 0.6498 \\ \\ \lambda_3^B &= -\lambda_1^B - \lambda_2^B = -0.7347 \\ \end{align*}\]

  4. Efek Interaksi: \[\begin{align*} \\ \lambda_{11}^{AB} &= \log(10) - \lambda - \lambda_1^A - \lambda_1^B \\ \\ &= 2.3026 - 2.4681 - 0.551 - 0.0425 = -0.759 \\ \\ \lambda_{12}^{AB} &= \log(25) - \lambda - \lambda_1^A - \lambda_2^B \\ \\ &= 3.2189 - 2.4681 - 0.551 - 0.6498 = -0.45 \\ \\ \lambda_{13}^{AB} &= \log(15) - \lambda - \lambda_1^A - \lambda_3^B \\ \\ &= 2.7081 - 2.4681 - 0.551 + 0.7347 = 0.4237 \\ \\ \lambda_{21}^{AB} &= \log(20) - \lambda - \lambda_2^A - \lambda_1^B \\ \\ &= 2.9957 - 2.4681 + 0.551 - 0.0425 = 1.0361 \\ \\ \lambda_{22}^{AB} &= \log(18) - \lambda - \lambda_2^A - \lambda_2^B \\ \\ &= 2.8904 - 2.4681 + 0.551 - 0.6494 = 0.3239 \\ \\ \lambda_{23}^{AB} &= \log(5) - \lambda - \lambda_2^A - \lambda_3^B \\ \\ &= 0.6931 - 2.4681 + 0.551 + 0.7347 = -0.4893 \\ \end{align*}\]

Hitung Odds Ratio dan Interval Kepercayaan

  1. Odds Ratio*: \[ \text{OR} = \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} = \frac{10 \cdot 18}{20 \cdot 25} = \frac{180}{500} = 0.36 \] \[ \log(\text{OR}) = \log\left( \frac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}} \right) = \log\left( \frac{180}{500} \right) = \log(0.36) = -1.0217 \]
  2. Hitung Standard Error:

\[ SE = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}}\\ \\ = \sqrt{\frac{1}{10} + \frac{1}{18} + \frac{1}{20} + \frac{1}{25}}\\ \\ = \sqrt{0.1 + 0.0556 + 0.05 + 0.04}\\ \\ = \sqrt{0.2456}\\ \\ = 0.4956 \\ \]

  1. Hitung Interval Kepercayaan 95%: \[\begin{align*} \\ \log(OR) \pm 1.96 \times SE \\ &= -1.0217 \pm 1.96 \times 0.4956\\ &= -1.0217 \pm 0.9714 = (-1.9931, -0.0503)\\ &= \text{CI}_{95\%} = (e^{-1.9931}, e^{-0.0503})\\ &= (0.1363, 0.9509) \\ \end{align*}\]

Jadi, OR=0.36 dengan taraf signifikansi 95% memiliki interval kepercayaan dengan rentang 0.1363 - 0.9509

Fitting Model Log-Linear di R

# Membuat data frame dari tabel
tabel2x3 <- matrix(c(10, 25, 15, 20, 18, 5), nrow = 2, byrow = TRUE)
colnames(tabel2x3) <- c("Rendah", "Sedang", "Tinggi")
rownames(tabel2x3) <- c("Pekerja", "Pengangguran")
tabel2x3
##              Rendah Sedang Tinggi
## Pekerja          10     25     15
## Pengangguran     20     18      5
# Ubah menjadi data.frame untuk glm
data2x3 <- as.data.frame(as.table(tabel2x3))
colnames(data2x3) <- c("Status", "TingkatKepuasan", "Freq")
data2x3
# Model log-linear tanpa interaksi (asumsi independen)
fit_no_inter23 <- glm(Freq ~ Status + TingkatKepuasan, family = poisson, data = data2x3)
summary(fit_no_inter23)
## 
## Call:
## glm(formula = Freq ~ Status + TingkatKepuasan, family = poisson, 
##     data = data2x3)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             2.7806     0.2064  13.475   <2e-16 ***
## StatusPengangguran     -0.1508     0.2080  -0.725    0.468    
## TingkatKepuasanSedang   0.3600     0.2379   1.513    0.130    
## TingkatKepuasanTinggi  -0.4055     0.2887  -1.405    0.160    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.4178  on 5  degrees of freedom
## Residual deviance:  9.2477  on 2  degrees of freedom
## AIC: 44.085
## 
## Number of Fisher Scoring iterations: 4
# Model log-linear dengan interaksi (untuk cek asosiasi)
fit_inter23 <- glm(Freq ~ Status * TingkatKepuasan, family = poisson, data = data2x3)
summary(fit_inter23)
## 
## Call:
## glm(formula = Freq ~ Status * TingkatKepuasan, family = poisson, 
##     data = data2x3)
## 
## Coefficients:
##                                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                2.3026     0.3162   7.281  3.3e-13
## StatusPengangguran                         0.6931     0.3873   1.790  0.07350
## TingkatKepuasanSedang                      0.9163     0.3742   2.449  0.01433
## TingkatKepuasanTinggi                      0.4055     0.4082   0.993  0.32062
## StatusPengangguran:TingkatKepuasanSedang  -1.0217     0.4955  -2.062  0.03924
## StatusPengangguran:TingkatKepuasanTinggi  -1.7918     0.6455  -2.776  0.00551
##                                             
## (Intercept)                              ***
## StatusPengangguran                       .  
## TingkatKepuasanSedang                    *  
## TingkatKepuasanTinggi                       
## StatusPengangguran:TingkatKepuasanSedang *  
## StatusPengangguran:TingkatKepuasanTinggi ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  1.8418e+01  on 5  degrees of freedom
## Residual deviance: -2.8866e-15  on 0  degrees of freedom
## AIC: 38.838
## 
## Number of Fisher Scoring iterations: 3

Interpretasi: Model tanpa interaksi: - Jika deviance tidak signifikan, maka Status Pekerjaan dan Tingkat Kepuasan Hidup independen.

  • Intercept: log frekuensi pada kategori referensi (Pekerja, Rendah)
  • Koefisien StatusPekerjaanPengangguran: perbedaan log-frekuensi antara Pengangguran vs Pekerja (pada tingkat kepuasan rendah)
  • Koefisien TingkatKepuasan: perbedaan log-frekuensi tingkat Kepuasan (Sedang/Tiggi) terhadap Rendah (pada status pekerjaan : pekerja)

Model dengan interaksi:

  • Jika koefisien interaksi (Status Pekerjaan:Tingkat Kepuasan) signifikan, berarti ada hubungan/asosiasi antara Status Pekerjaan dan Tingkat Kepuasan. Artinya distribusi Tingkat Kepuasan berbeda antara Pekerja dan Pengangguran.

Contoh interpretasi hasil (misal): - Jika koefisien StatusPekerjaanPengangguran negatif: proporsi Pengangguran pada kategori referensi lebih kecil dibanding Pekerja. - Jika koefisien TingkatKepuasan_Sedang positif: kemungkinan seseorang dengan tingkat kepuasan Sedang lebih tinggi daripada Rendah (pada Pekerja). - Jika model interaksi signifikan, pola distribusi Tingkat Kepuasan berbeda antara Pekerja dan Pengangguran.

Model Log Linear Tiga Arah

Pada bagian sebelumnya, telah dijelaskan bahwa salah satu tujuan utama dari penggunaan model log-linear adalah untuk memperkirakan parameter-parameter yang menggambarkan hubungan antar variabel kategorik.

Dalam materi kali ini, kita akan mempelajari model log-linear yang lebih kompleks, yaitu yang digunakan untuk tabel kontingensi tiga arah. Model ini melibatkan tiga variabel kategorik, sehingga kemungkinan munculnya interaksi antar variabel menjadi lebih banyak. Dalam hal ini, interaksi paling kompleks yang dapat dianalisis adalah interaksi tiga arah, yaitu interaksi yang melibatkan ketiga variabel sekaligus.

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

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

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

  3. 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}\)

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

Contoh Kasus

Deskripsi Data

Sebuah tim peneliti psikologi dari universitas ternama ingin memahami bagaimana faktor-faktor seperti jenis kelamin dan frekuensi olahraga berkaitan dengan tingkat stres pada individu. Stres merupakan salah satu isu kesehatan mental yang umum, terutama di kalangan mahasiswa dan pekerja muda, sehingga penting untuk memahami faktor-faktor yang mungkin berkontribusi terhadap tingginya atau rendahnya tingkat stres seseorang.

Peneliti mengklasifikasikan responden ke dalam tiga dimensi utama: Tingkat Stres: - Tinggi - Sedang - Rendah

Jenis Kelamin: - Laki-laki - Perempuan

Frekuensi Olahraga: - Jarang (kurang dari 1x per minggu) - Sedang (1-2x per minggu) - Rutin (lebih dari 3x per minggu)

Survei dilakukan kepada 462 responden dari berbagai latar belakang. Data dikumpulkan dan disusun dalam format tiga dimensi, dan berikut adalah salah satu cara menyajikannya, yaitu dalam bentuk tabel silang per jenis kelamin.

Tabel data:

Tingkat Stres Jenis Kelamin Jarang Rutin Total
Tinggi Laki-laki 52 20 72
Perempuan 60 16 76
Total 112 36 148
Sedang Laki-laki 30 38 68
Perempuan 36 36 72
Total 66 74 140
Rendah Laki-laki 18 62 80
Perempuan 14 80 94
Total 32 142 174

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.stres <- factor(rep(c("1tinggi", "2sedang", "3rendah"), each = 4))
x.sex  <- factor(rep(c("1L", "2P"), each = 2, times = 3))
y.freq  <- factor(rep(c("1jarang", "2rutin"), times = 6))
counts <- c(52, 20, 60, 16, 30, 38, 36, 36, 18, 62, 14, 80)

data3 <- data.frame(
  TingkatStres              = z.stres,
  JenisKelamin               = x.sex,
  FrekuensiOlahraga          = y.freq,
  Frekuensi                  = counts

)
data3

Membentuk Tabel Kontingensi 3 Arah

table3d <- xtabs(Frekuensi ~ TingkatStres + JenisKelamin + FrekuensiOlahraga, data = data3)
ftable(table3d)
##                           FrekuensiOlahraga 1jarang 2rutin
## TingkatStres JenisKelamin                                 
## 1tinggi      1L                                  52     20
##              2P                                  60     16
## 2sedang      1L                                  30     38
##              2P                                  36     36
## 3rendah      1L                                  18     62
##              2P                                  14     80

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 = "2P")
y.freq  <- relevel(y.freq, ref = "2rutin")
z.stres <- relevel(z.stres, ref = "3rendah")

Model Saturated memasukkan semua interaksi hingga tiga arah (tiga kategori): \[ \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.freq + z.stres +
             x.sex*y.freq + x.sex*z.stres + y.freq*z.stres +
             x.sex*y.freq*z.stres,
             family = poisson(link = "log"))
summary(model_saturated)
## 
## Call:
## glm(formula = counts ~ x.sex + y.freq + z.stres + x.sex * y.freq + 
##     x.sex * z.stres + y.freq * z.stres + x.sex * y.freq * z.stres, 
##     family = poisson(link = "log"))
## 
## Coefficients:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                            4.3820     0.1118  39.194  < 2e-16 ***
## x.sex1L                               -0.2549     0.1692  -1.506    0.132    
## y.freq1jarang                         -1.7430     0.2897  -6.016 1.78e-09 ***
## z.stres1tinggi                        -1.6094     0.2739  -5.877 4.18e-09 ***
## z.stres2sedang                        -0.7985     0.2007  -3.979 6.93e-05 ***
## x.sex1L:y.freq1jarang                  0.5062     0.3945   1.283    0.199    
## x.sex1L:z.stres1tinggi                 0.4780     0.3757   1.272    0.203    
## x.sex1L:z.stres2sedang                 0.3090     0.2876   1.074    0.283    
## y.freq1jarang:z.stres1tinggi           3.0647     0.4039   7.589 3.23e-14 ***
## y.freq1jarang:z.stres2sedang           1.7430     0.3735   4.667 3.06e-06 ***
## x.sex1L:y.freq1jarang:z.stres1tinggi  -0.8725     0.5514  -1.582    0.114    
## x.sex1L:y.freq1jarang:z.stres2sedang  -0.7426     0.5204  -1.427    0.154    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1.2499e+02  on 11  degrees of freedom
## Residual deviance: 2.5757e-14  on  0  degrees of freedom
## AIC: 88.183
## 
## Number of Fisher Scoring iterations: 3
exp(model_saturated$coefficients)
##                          (Intercept)                              x.sex1L 
##                           80.0000000                            0.7750000 
##                        y.freq1jarang                       z.stres1tinggi 
##                            0.1750000                            0.2000000 
##                       z.stres2sedang                x.sex1L:y.freq1jarang 
##                            0.4500000                            1.6589862 
##               x.sex1L:z.stres1tinggi               x.sex1L:z.stres2sedang 
##                            1.6129032                            1.3620072 
##         y.freq1jarang:z.stres1tinggi         y.freq1jarang:z.stres2sedang 
##                           21.4285714                            5.7142857 
## x.sex1L:y.freq1jarang:z.stres1tinggi x.sex1L:y.freq1jarang:z.stres2sedang 
##                            0.4179259                            0.4758772

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), frekuensi olahraga (y.freq), dan tingkat stress (z.stres) terhadap frekuensi responden. Dengan hasil estimasi koefisien nya sebagai berikut:

Parameter Estimate Std. Error z value Pr(>)
(Intercept) 4.3820 0.1118 39.194 <2e-16 ***
x.sex1L -0.2549 0.1692 -1.506 0.132
y.freq1jarang -1.7430 0.2897 -6.016 1.78e-09 ***
z.stres1tinggi -1.6094 0.2739 -5.877 4.18e-09 ***
z.stres2sedang -0.7985 0.2007 -3.979 6.93e-05 ***
x.sex1L:y.freq1jarang 0.5062 0.3945 1.283 0.199
x.sex1L:z.stres1tinggi 0.4780 0.3757 1.272 0.203
x.sex1L:z.stres2sedang 0.3090 0.2876 1.074 0.283
y.freq1jarang:z.stres1tinggi 3.0647 0.4309 7.108 3.23e-14 ***
y.freq1jarang:z.stres2sedang 1.7430 0.3735 4.667 3.06e-06 ***
x.sex1L:y.freq1jarang:z.stres1tinggi -0.8725 0.5514 -1.582 0.114
x.sex1L:y.freq1jarang:z.stres2sedang -0.7426 0.5204 -1.427 0.154

Interpretasi Koefisien:

  • (Intercept):
    Rata-rata log jumlah kasus untuk kategori referensi (Perempuan, Frekuensi Sering, Stres Rendah) adalah 4.382. Ini menunjukkan baseline expected count sebesar exp(4.382) ≈ 80.

  • x.sex1L (Laki-laki):
    Laki-laki memiliki expected count sekitar exp(-0.2549) ≈ 0.775 kali perempuan, tetapi perbedaannya tidak signifikan (p = 0.132).

  • y.freq1jarang (Frekuensi olahraga jarang):
    Responden dengan frekuensi olahraga yang jarang memiliki expected count sekitar exp(-1.7430) ≈ 0.175 dari mereka yang sering (signifikan, p < 0.001).

  • z.stres1tinggi (Stres tinggi):
    Responden dengan tingkat stres tinggi memiliki expected count sekitar exp(-1.6094) ≈ 0.20 dari mereka yang stres rendah (signifikan, p < 0.001).

  • z.stres2sedang (Stres sedang):
    Responden dengan stres sedang memiliki expected count sekitar exp(-0.7985) ≈ 0.45 dari yang stres rendah (signifikan, p < 0.001).

  • x.sex1L:y.freq1jarang:
    Interaksi ini menunjukkan bahwa efek frekuensi olahraga jarang pada laki-laki sedikit lebih tinggi (exp(0.5062) ≈ 1.66), tetapi tidak signifikan (p = 0.199).

  • x.sex1L:z.stres1tinggi:
    Interaksi ini menunjukkan bahwa stres tinggi pada laki-laki meningkatkan expected count (exp(0.4780) ≈ 1.61), tapi tidak signifikan (p = 0.203).

  • x.sex1L:z.stres2sedang:
    Interaksi stres sedang pada laki-laki (exp(0.3090) ≈ 1.36) juga tidak signifikan (p = 0.283).

  • y.freq1jarang:z.stres1tinggi:
    Interaksi antara frekuensi olahraga jarang dan stres tinggi meningkatkan expected count sebesar exp(3.0647) ≈ 21.4, sangat signifikan (p < 0.001).

  • y.freq1jarang:z.stres2sedang:
    Interaksi antara olahraga jarang dan stres sedang menghasilkan expected count exp(1.7430) ≈ 5.7 kali lipat, signifikan (p < 0.001).

  • x.sex1L:y.freq1jarang:z.stres1tinggi:
    Interaksi tiga arah ini menurunkan expected count exp(-0.8725) ≈ 0.42, tetapi tidak signifikan (p = 0.114).

  • x.sex1L:y.freq1jarang:z.stres2sedang:
    Interaksi tiga arah ini juga menurunkan expected count exp(-0.7426) ≈ 0.48, tidak signifikan (p = 0.154).

Goodness-of-Fit

  • Residual deviance ≈ 0:
    Menandakan bahwa model saturated sangat cocok dengan data (seluruh variasi data dijelaskan model).

  • AIC = 88.18:
    Dapat digunakan untuk perbandingan dengan model yang lebih sederhana. Semakin kecil AIC, semakin baik model.

Model ini sangat fit terhadap data (residual deviance mendekati nol), tetapi tidak semua parameter signifikan.

Kesimpulan:

  • Model ini sangat fit terhadap data (residual deviance mendekati nol), tetapi tidak semua parameter signifikan.

  • Efek utama yang paling signifikan adalah:

  • Frekuensi olahraga jarang (expected count turun drastis dibandingkan yang sering)

  • Tingkat stres tinggi dan sedang (expected count turun signifikan dibandingkan stres rendah)

  • Interaksi antara olahraga jarang dan stres tinggi/sedang meningkatkan expected count secara signifikan

  • Tidak ditemukan bukti kuat untuk interaksi 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.freq + z.stres +
              x.sex*y.freq + x.sex*z.stres + y.freq*z.stres,
              family = poisson(link = "log"))
summary(model_homogenous)
## 
## Call:
## glm(formula = counts ~ x.sex + y.freq + z.stres + x.sex * y.freq + 
##     x.sex * z.stres + y.freq * z.stres, family = poisson(link = "log"))
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   4.33523    0.11098  39.064  < 2e-16 ***
## x.sex1L                      -0.15079    0.15706  -0.960 0.337017    
## y.freq1jarang                -1.46409    0.21804  -6.715 1.88e-11 ***
## z.stres1tinggi               -1.43945    0.22433  -6.417 1.39e-10 ***
## z.stres2sedang               -0.70931    0.18283  -3.880 0.000105 ***
## x.sex1L:y.freq1jarang        -0.05711    0.21327  -0.268 0.788866    
## x.sex1L:z.stres1tinggi        0.13993    0.25524   0.548 0.583531    
## x.sex1L:z.stres2sedang        0.12054    0.23561   0.512 0.608923    
## y.freq1jarang:z.stres1tinggi  2.62706    0.27403   9.587  < 2e-16 ***
## y.freq1jarang:z.stres2sedang  1.37740    0.25888   5.321 1.03e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 124.9856  on 11  degrees of freedom
## Residual deviance:   2.9849  on  2  degrees of freedom
## AIC: 87.168
## 
## Number of Fisher Scoring iterations: 4

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)

  1. Hitung Selisih Deviance
# Deviance antar model
Deviance.model <- model_homogenous$deviance - model_saturated$deviance
Deviance.model
## [1] 2.984885
  1. 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
  1. Chi-Square Tabel (\(\alpha\)=0.05)
chi.tabel <- qchisq(1 - 0.05, df = derajat.bebas)
chi.tabel
## [1] 5.991465
  1. 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, frekuensi olahraga, dan tingkat stress.

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 = 2.9849 - 0.00 = 2.9849

    • \(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 2.9849 < 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, frekuensi olahraga, dan tingkat stress.

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.freq + z.stres +
              x.sex*y.freq + x.sex*z.stres,
              family = poisson(link = "log"))
summary(model_conditional_X)
## 
## Call:
## glm(formula = counts ~ x.sex + y.freq + z.stres + x.sex * y.freq + 
##     x.sex * z.stres, family = poisson(link = "log"))
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             3.937e+00  1.187e-01  33.178   <2e-16 ***
## x.sex1L                -1.613e-01  1.743e-01  -0.925   0.3548    
## y.freq1jarang          -1.823e-01  1.291e-01  -1.412   0.1579    
## z.stres1tinggi         -2.126e-01  1.543e-01  -1.378   0.1682    
## z.stres2sedang         -2.666e-01  1.566e-01  -1.702   0.0887 .  
## x.sex1L:y.freq1jarang  -3.517e-15  1.871e-01   0.000   1.0000    
## x.sex1L:z.stres1tinggi  1.072e-01  2.240e-01   0.479   0.6323    
## x.sex1L:z.stres2sedang  1.041e-01  2.274e-01   0.458   0.6471    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 124.99  on 11  degrees of freedom
## Residual deviance: 115.77  on  4  degrees of freedom
## AIC: 195.95
## 
## Number of Fisher Scoring iterations: 5

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 frekuensi olahraga (Y) dan tingkat stress (Z)) -

    • H0: \(\lambda_{jk}^{YZ}\neq0\) (Ada interaksi antara frekuensi olahraga (Y) dan tingkat stress (Z))

  • Tingkat Signifikansi

    \(\alpha=5\%\)

  • Statistik Uji

    • \(\Delta\)Deviance = Deviance model conditional on X - Deviance model homogenous = 115.77 - 2.9849 = 112.7855

    • \(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 112.7855 < 5.991, maka tolak H0

  • Interpretasi

    Dengan taraf nyata 5%, sudah cukup bukti untuk menolak H0 atau dapat dikatakan bahwa ada interaksi antara frekuensi olahraga dan tingkat stress. Dengan kata lain, model yang terbentuk adalah model dengan 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] 112.7855
  1. Hitung Derajat Bebas
# Chi Square tabel dengan alpha = 0.05
derajat.bebas <- (4 - 2)
derajat.bebas
## [1] 2
  1. Chi-Square Tabel
chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 5.991465
  1. Keputusan Uji
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan
## [1] "Tolak"

Interpretasi Karena nilai Deviance.model = 112.7855 lebih kecil dari nilai kritis chi-square tabel = 5.99 (dengan df = 2, alpha = 0.05), maka keputusan uji adalah “Tolak H0”.

Pada taraf nyata 5%, sudah cukup bukti untuk menolak H0, atau dengan kata lain ada interaksi antara frekuesi olahraga (Y) dan tingkat stress (Z). Model yang terbentuk adalah model homogenous, sehingga interaksi Y*Z 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.freq + z.stres +
              x.sex*y.freq + y.freq*z.stres,
              family = poisson(link = "log"))
summary(model_conditional_Y)
## 
## Call:
## glm(formula = counts ~ x.sex + y.freq + z.stres + x.sex * y.freq + 
##     y.freq * z.stres, family = poisson(link = "log"))
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   4.309e+00  1.032e-01  41.757  < 2e-16 ***
## x.sex1L                      -9.531e-02  1.261e-01  -0.756     0.45    
## y.freq1jarang                -1.490e+00  2.150e-01  -6.930 4.20e-12 ***
## z.stres1tinggi               -1.372e+00  1.866e-01  -7.354 1.92e-13 ***
## z.stres2sedang               -6.518e-01  1.434e-01  -4.546 5.47e-06 ***
## x.sex1L:y.freq1jarang         2.851e-15  1.871e-01   0.000     1.00    
## y.freq1jarang:z.stres1tinggi  2.625e+00  2.739e-01   9.585  < 2e-16 ***
## y.freq1jarang:z.stres2sedang  1.376e+00  2.588e-01   5.316 1.06e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 124.9856  on 11  degrees of freedom
## Residual deviance:   3.3587  on  4  degrees of freedom
## AIC: 83.541
## 
## 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 tingkat stress (Z))

    • H0: \(\lambda_{ik}^{XZ}\neq0\) (Ada interaksi antara jenis kelamin

      1. dan stress (Z))
  • Tingkat Signifikansi

    \(\alpha=5\%\)

  • Statistik Uji

    • \(\Delta\)Deviance = Deviance model conditional on Y - Deviance model homogenous = 3.3597 - 2.9849 = 0.3738

    • \(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 0.3738 < 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 tingkat stress. 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] 0.3737977
  1. Hitung Derajat Bebas
# Chi Square tabel dengan alpha = 0.05
derajat.bebas <- (4 - 2)
derajat.bebas
## [1] 2
  1. Chi-Square Tabel
chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 5.991465
  1. Keputusan Uji
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan
## [1] "Terima"

Interpretasi Karena nilai Deviance.model = 0.3738 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 tingkat stress (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.freq + z.stres +
              x.sex*z.stres + y.freq*z.stres,
              family = poisson(link = "log"))
summary(model_conditional_Z)
## 
## Call:
## glm(formula = counts ~ x.sex + y.freq + z.stres + x.sex * z.stres + 
##     y.freq * z.stres, family = poisson(link = "log"))
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    4.3401     0.1092  39.730  < 2e-16 ***
## x.sex1L                       -0.1613     0.1521  -1.060    0.289    
## y.freq1jarang                 -1.4901     0.1957  -7.615 2.64e-14 ***
## z.stres1tinggi                -1.4230     0.2147  -6.627 3.43e-11 ***
## z.stres2sedang                -0.7010     0.1794  -3.907 9.35e-05 ***
## x.sex1L:z.stres1tinggi         0.1072     0.2240   0.479    0.632    
## x.sex1L:z.stres2sedang         0.1041     0.2274   0.458    0.647    
## y.freq1jarang:z.stres1tinggi   2.6251     0.2739   9.585  < 2e-16 ***
## y.freq1jarang:z.stres2sedang   1.3757     0.2588   5.316 1.06e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 124.9856  on 11  degrees of freedom
## Residual deviance:   3.0566  on  3  degrees of freedom
## AIC: 85.239
## 
## 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 frekuensi olahraga (Y)) -

    • H0: \(\lambda_{ij}^{XY}\neq0\) (Ada interaksi antara jenis kelamin

      1. dan frekuesni olahraga (Y))
  • Tingkat Signifikansi

    \(\alpha=5\%\)

  • Statistik Uji

    • \(\Delta\)Deviance = Deviance model conditional on Y - Deviance model homogenous = 3.0566 - 2.9849 = 0.072

    • \(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 0.072 < 3.841, maka tolak H0

  • Interpretasi

    Dengan taraf nyata 5%, menolak H0 atau dapat dikatakan bahwa ada interaksi antara jenis kelamin dan frekuensi olahraga. 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] 0.07173518
  1. Hitung Derajat Bebas
# Chi Square tabel dengan alpha = 0.05
derajat.bebas <- (3 - 2)
derajat.bebas
## [1] 1
  1. Chi-Square Tabel
chi.tabel <- qchisq((1 - 0.05), df = derajat.bebas)
chi.tabel
## [1] 3.841459
  1. Keputusan Uji
Keputusan <- ifelse(Deviance.model <= chi.tabel, "Terima", "Tolak")
Keputusan
## [1] "Terima"

Interpretasi Karena nilai Deviance.model = 0.072 lebih kecil dari nilai kritis chi-square tabel = 3.84 (dengan df = 1, alpha = 0.05), maka keputusan uji adalah “Terima H0”.

Pada taraf nyata 5%, belum cukup bukti untuk menolak H0, atau dengan kata lain tidak ada interaksi antara jenis kelamin (X) dan frekuesni olahraga (Y). Model yang terbentuk cukup hanya sampai dua interaksi dengan Z (conditional on Z), sehingga interaksi XY tidak signifikan secara statistik.

Pemilihan Model Terbaik

Ringkasan Model Log Linear

Model Parameter Deviance Jumlah Parameter df AIC
Saturated λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢⱼˣʸ + λᵢₖˣᶻ + λⱼₖʸᶻ + λᵢⱼₖˣʸᶻ 0.000 12 0 88.18
Homogenous λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢⱼˣʸ + λᵢₖˣᶻ + λⱼₖʸᶻ 2.9849 10 2 87.168
Conditional on X λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢⱼˣʸ + λᵢₖˣᶻ 115.77 8 4 195.95
Conditional on Y λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢⱼˣʸ + λⱼₖʸᶻ 3.3597 8 4 83.541
Conditional on Z λ + λᵢˣ + λⱼʸ + λₖᶻ + λᵢₖˣᶻ + λⱼₖʸᶻ 3.0566 9 3 85.239

Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah

Interaksi Pengujian \(\Delta\)Deviance \(\Delta\)df Chi-Square Tabel Keputusan Keterangan
XYZ Saturated vs Homogenous 2.9849 2 5.991 Terima H0 tidak ada interaksi
YZ Conditional on X vs Homogenous 112.7855 2 5.991 Tolak H0 ada interaksi
XZ Conditional on Y vs Homogenous 0.3738 2 5.991 Terima H0 tidak ada interaksi
XY Conditional on Z vs Homogenous 0.072 1 3.841 Terima H0 tidak ada interaksi

Kesimpulan Pemilihan Model Terbaik

Dari hasil di atas diketahui bahwa asosiasi yang nyata hanya terdapat antara frekuensi olahraga dan tingkat stress. Sehingga, model terbaik adalah:

\[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{YZ}_{jk} \]

Model terbaik adalah model log-linear tanpa interaksi tiga arah dan hanya memuat interaksi dua arah antara frekuesni olahraga (Y) dan tingkat stress (Z).

Model Terbaik

Model terbaik dipilih berdasarkan pengujian interaksi yang signifikan, yaitu hanya interaksi dua arah antara frekuensi olahraga (Y) dan tingkat stress (Z):

\[ \log(\mu_{ijk}) = \lambda + \lambda^{X}_i + \lambda^{Y}_j + \lambda^{Z}_k + \lambda^{YZ}_{jk} \]

# Model Terbaik
bestmodel <- glm(counts ~ x.sex + y.freq + z.stres +
                 y.freq*z.stres,
                 family = poisson(link = "log"))
summary(bestmodel)
## 
## Call:
## glm(formula = counts ~ x.sex + y.freq + z.stres + y.freq * z.stres, 
##     family = poisson(link = "log"))
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   4.30920    0.09492  45.398  < 2e-16 ***
## x.sex1L                      -0.09531    0.09315  -1.023    0.306    
## y.freq1jarang                -1.49009    0.19568  -7.615 2.64e-14 ***
## z.stres1tinggi               -1.37231    0.18660  -7.354 1.92e-13 ***
## z.stres2sedang               -0.65176    0.14337  -4.546 5.47e-06 ***
## y.freq1jarang:z.stres1tinggi  2.62507    0.27386   9.585  < 2e-16 ***
## y.freq1jarang:z.stres2sedang  1.37568    0.25876   5.316 1.06e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 124.9856  on 11  degrees of freedom
## Residual deviance:   3.3587  on  5  degrees of freedom
## AIC: 81.541
## 
## Number of Fisher Scoring iterations: 4
Model AIC
Saturated 88.18
Homogenous 87.168
Conditional on X 195.95
Conditional on Y 83.541
Conditional on Z 85.239
Model Terbaik 81.541

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.309 74.3809
x.sex1M -0.0953 0.9091
y.freq1jarang -1.49 0.225
z.stress1tinggi -1.372 0.2535
z.stress2sedang -0.652 0.5211
z.stress1tinggi:y.freq1jarang 2.625 13.806

Interpretasi Koefisien Model Terbaik

  • \(\exp(\lambda_{\text{1M}}^X) = \exp(-0.0953) = 0.909 \rightarrow\) nilai odds

    Tanpa memperhatikan frekuensi olahraga dan tingkat stres, peluang seseorang berjenis kelamin laki-laki adalah 0,909 kali dibandingkan perempuan.
    Atau, peluang seseorang perempuan adalah \(1 / 0.909 = 1.10\) kali dibandingkan laki-laki.

  • \(\exp(\lambda_{\text{1jarang}}^Y) = \exp(-1.49) = 0.225 \rightarrow\) nilai odds

    Tanpa memperhatikan jenis kelamin dan tingkat stres, peluang seseorang yang berolahraga jarang adalah 0.225 kali dibandingkan yang sering.
    Atau, peluang seseorang yang sering berolahraga adalah \(1 / 0.225 = 4.44\) kali dibandingkan yang jarang.

  • \(\exp(\lambda_{\text{1tinggi}}^Z) = \exp(-1.372) = 0.253 \rightarrow\) nilai odds

    Tanpa memperhatikan jenis kelamin dan frekuensi olahraga, peluang seseorang dengan tingkat stres tinggi adalah 0.253 kali dibandingkan yang tingkat stresnya rendah.
    Atau, peluang seseorang dengan stres rendah adalah \(1 / 0.253 = 3.95\) kali dibandingkan dengan yang stres tinggi.

  • \(\exp(\lambda_{\text{2sedang}}^Z) = \exp(-0.652) = 0.521 \rightarrow\) nilai odds

    Tanpa memperhatikan jenis kelamin dan frekuensi olahraga, peluang seseorang dengan tingkat stres sedang adalah 0.521 kali dibandingkan yang rendah.
    Atau, peluang seseorang dengan stres rendah adalah \(1 / 0.521 = 1.92\) kali dibandingkan yang stres sedang.

  • \(\exp(\lambda_{\text{1tinggi,1jarang}}^{Z,Y}) = \exp(2.625) = 13.806 \rightarrow\) nilai odds ratio

    Tanpa memperhatikan jenis kelamin, odds seseorang yang memiliki stres tinggi dan olahraga jarang adalah 13.806 kali lebih tinggi dibandingkan dengan orang yang tidak memiliki kedua kondisi tersebut.

Nilai Dugaan Model Terbaik

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

\[\begin{align*} \\ \hat{\mu}_{111} &= \exp\left( \lambda + \lambda^{x}_{1m} + \lambda^{y}_{1jarang} + \lambda^{z}_{1tinggi} + \lambda^{yz}_{1jarang,1tinggi} \right) \\ \\ &= \exp\left( 4.309 - 0.0953 - 1.49 - 1.372 + 2.625 \right) \\ \\ &= \exp(3.9767) = 53.333 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{112} &= \exp\left( \lambda + \lambda^{x}_{1m} + \lambda^{y}_{1jarang} + \lambda^{z}_{2sedang} + \lambda^{yz}_{1jarang,2sedang} \right) \\ \\ &= \exp\left( 4.309 - 0.0953 - 1.49 - 0.652 + 1.376 \right) \\ \\ &= \exp(3.4477) = 31.429 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{113} &= \exp\left( \lambda + \lambda^{x}_{1m} + \lambda^{y}_{1jarang} + \lambda^{z}_{3rendah} + \lambda^{yz}_{1jarang,3rendah} \right) \\ \\ &= \exp\left( 4.309 - 0.0953 - 1.49 + 0 + 0 \right) \\ \\ &= \exp(2.7237) = 15.2381 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{121} &= \exp\left( \lambda + \lambda^{x}_{1m} + \lambda^{y}_{2rutin} + \lambda^{z}_{1tinggi} + \lambda^{yz}_{2rutin,1tinggi} \right) \\ \\ &= \exp\left( 4.309 - 0.0953 + 0 - 1.372 + 0 \right) \\ \\ &= \exp(2.8417) = 17.143 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{122} &= \exp\left( \lambda + \lambda^{x}_{1m} + \lambda^{y}_{2rutin} + \lambda^{z}_{2sedang} + \lambda^{yz}_{2rutin,2sedang} \right) \\ \\ &= \exp\left( 4.309 - 0.0953 + 0 - 0.652 + 0 \right) \\ \\ &= \exp(3.5617) = 35.238 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{123} &= \exp\left( \lambda + \lambda^{x}_{1m} + \lambda^{y}_{2rutin} + \lambda^{z}_{3rendah} + \lambda^{yz}_{2rutin,3rendah} \right) \\ \\ &= \exp\left( 4.309 - 0.0953 + 0 + 0 + 0 \right) \\ \\ &= \exp(4.2137) = 67.619 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{211} &= \exp\left( \lambda + \lambda^{x}_{2f} + \lambda^{y}_{1jarang} + \lambda^{z}_{1tinggi} + \lambda^{yz}_{1jarang,1tinggi} \right) \\ \\ &= \exp\left( 4.309 + 0 - 1.49 - 1.372 + 2.625 \right) \\ \\ &= \exp(4.072) = 58.667 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{212} &= \exp\left( \lambda + \lambda^{x}_{2f} + \lambda^{y}_{1jarang} + \lambda^{z}_{2sedang} + \lambda^{yz}_{1jarang,2sedang} \right) \\ \\ &= \exp\left( 4.309 + 0 - 1.49 - 0.652 + 1.376 \right) \\ \\ &= \exp(3.543) = 34.571 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{213} &= \exp\left( \lambda + \lambda^{x}_{2f} + \lambda^{y}_{1jarang} + \lambda^{z}_{3rendah} + \lambda^{yz}_{1jarang,3rendah} \right) \\ \\ &= \exp\left( 4.309 + 0 - 1.49 + 0 + 0 \right) \\ \\ &= \exp(2.819) = 16.762 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{221} &= \exp\left( \lambda + \lambda^{x}_{2f} + \lambda^{y}_{2rutin} + \lambda^{z}_{1tinggi} + \lambda^{yz}_{2rutin,1tinggi} \right) \\ \\ &= \exp\left( 4.309 + 0 + 0 - 1.372 + 0 \right) \\ \\ &= \exp(2.937) = 18.857 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{222} &= \exp\left( \lambda + \lambda^{x}_{2f} + \lambda^{y}_{2rutin} + \lambda^{z}_{2sedang} + \lambda^{yz}_{2rutin,2sedang} \right) \\ \\ &= \exp\left( 4.309 + 0 + 0 - 0.652 + 0 \right) \\ \\ &= \exp(3.657) = 38.762 \\ \end{align*}\]

\[\begin{align*} \\ \hat{\mu}_{223} &= \exp\left( \lambda + \lambda^{x}_{2f} + \lambda^{y}_{2rutin} + \lambda^{z}_{3rendah} + \lambda^{yz}_{2rutin,3rendah} \right) \\ \\ &= \exp\left( 4.309 + 0 + 0 + 0 + 0 \right) \\ \\ &= \exp(4.309) = 74.381 \\ \end{align*}\]

Keterangan: Nilai \(\hat{\mu}_{ijk}\) akan sama apapun referensi dari kategori peubahnya yang digunakan.

Dengan menggunakan R:

# Fitted values dari model terbaik
data.frame(
  Stress = z.stres,
  sex = x.sex,
  Freq = y.freq,
  counts = counts,
  fitted = bestmodel$fitted.values
)

REFRENSI