KATA PENGANTAR
I. Tabel Kontingensi 2 x 2
II. Inferensi Tabel Kontingensi Dua Arah
III. Tabel Kontingensi Tiga Arah
IV. Inferensi Tabel Kontingensi Tiga Arah
- Pengujian Independensi
V. Generalized Linear Model (GLM)
VI. Inferensi Generalized Linear Model (GLM)
- Inferensi Parameter
VII. Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio
- Contoh Kasus
VIII. Pemilihan Model Regresi Logistik dan Evaluasi
IX. Distribusi Multinomial
X. Regresi Logistik Ordinal
XI. Log Linear Model
XII. Referensi

KATA PENGANTAR

Puji dan syukur penulis panjatkan ke hadirat Tuhan Yang Maha Esa atas segala rahmat dan karunia-Nya, sehingga penulis dapat menyusun dan menyelesaikan eBook ini sebagai bagian dari pembelajaran dalam mata kuliah Analisis Data Kategorik. Buku ini disusun dengan tujuan untuk memberikan pemahaman yang sistematis dan aplikatif mengenai berbagai metode analisis yang digunakan untuk menangani data kategorik, baik dalam konteks dua variabel maupun dalam struktur multivariat yang lebih kompleks.

Isi buku ini mencakup beragam topik penting, antara lain tabel kontingensi dua arah (2×2) dan perluasannya ke dalam bentuk tiga arah beserta inferensinya, model linear umum terumumkan (Generalized Linear Model), regresi logistik biner, multinomial, dan ordinal, strategi pemilihan model regresi logistik, distribusi multinomial, serta model log-linear. Seluruh materi disajikan dengan mengedepankan keseimbangan antara teori dan praktik, sehingga pembaca diharapkan tidak hanya memahami dasar-dasar konsep statistik yang digunakan, tetapi juga mampu menerapkannya dalam analisis data nyata.

Penyusunan eBook ini diharapkan dapat menjadi referensi yang bermanfaat bagi mahasiswa, dosen, maupun praktisi yang tertarik dalam analisis data kategorik, baik di bidang kesehatan, ilmu sosial, ekonomi, maupun sains terapan lainnya. Penulis menyadari bahwa buku ini masih jauh dari sempurna. Oleh karena itu, segala bentuk kritik dan saran yang membangun sangat penulis harapkan untuk penyempurnaan di masa yang akan datang. Akhir kata, semoga eBook ini dapat memberikan manfaat dan menjadi bagian dari upaya peningkatan literasi statistik di Indonesia.

I. Tabel Kontingensi 2 x 2

Tabel kontingensi 2 x 2 adalah bentuk sederhana dari tabel silang yang digunakan untuk menyajikan dan menganalisis hubungan antara dua variabel kategorik, masing-masing dengan dua kategori. Tabel ini menampilkan frekuensi kombinasi dari pasangan kategori kedua variabel, sehingga memungkinkan untuk menghitung peluang bersama, peluang marginal, dan peluang bersyarat. Melalui struktur yang ringkas, tabel kontingensi 2 x 2 membantu mengidentifikasi ada tidaknya keterkaitan atau ketergantungan antara dua variabel. Tabel ini juga menjadi dasar dalam berbagai analisis statistik inferensial, seperti uji independensi atau pengukuran asosiasi antara variabel.

Tabel kontingensi 2 x 2 memiliki struktur sebagai berikut:

	Kategori 1 (+)	Kategori 2 (-)	Total
Grup 1	\(n_{11}\)	\(n_{12}\)	\(n_{1.}\)
Grup 2	\(n_{21}\)	\(n_{22}\)	\(n_{2.}\)
Total	\(n_{.1}\)	\(n_{.2}\)	\(n\)

Contoh Kasus

Sebuah perusahaan ingin mengetahui apakah terdapat hubungan antara mengikuti pelatihan kerja dan status kelulusan ujian sertifikasi. Dari 120 karyawan yang mengikuti ujian, diperoleh data sebagai berikut:

70 orang mengikuti pelatihan, sisanya tidak.
Dari yang mengikuti pelatihan, 56 orang lulus ujian, dan 14 tidak lulus.
Dari yang tidak mengikuti pelatihan, 18 orang lulus dan 32 tidak lulus.

Data tersebut dapat disajikan dalam tabel kontingensi 2x2 berikut:

	Lulus	Tidak Lulus	Total
Ikut Pelatihan	56	14	70
Tidak Ikut Pelatihan	18	32	50
Total	74	46	120

Notasi

A : Mengikuti pelatihan
A’ : Tidak mengikuti pelatihan
B : Lulus ujian
B’ : Tidak lulus ujian

Distribusi Peluang

- Peluang Bersama (Joint Probability)

Probabilitas terjadinya dua kejadian secara bersamaan, dihitung dengan membagi frekuensi sel tertentu dengan total keseluruhan data.

Rumus umum:

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

Perhitungan Manual:

P(Ikut Pelatihan \(\cap\) Lulus) \(= \frac{56}{120} = 0.467\)
P(Ikut Pelatihan \(\cap\) Tidak Lulus) \(= \frac{14}{120} = 0.117\)
P(Tidak Ikut Pelatihan \(\cap\) Lulus) \(= \frac{18}{120} = 0.15\)
P(Tidak Ikut Pelatihan \(\cap\) Tidak Lulus) \(= \frac{32}{120} = 0.267\)

Perhitungan dengan R:

knitr::opts_chunk$set(echo = TRUE)
library(tibble)

# Data frekuensi
total <- 120
a <- 56   # Ikut pelatihan dan lulus
b <- 14   # Ikut pelatihan dan tidak lulus
c <- 18   # Tidak ikut pelatihan dan lulus
d <- 32   # Tidak ikut pelatihan dan tidak lulus

# Peluang bersama
p_AB   <- a / total
p_ABp  <- b / total
p_ApB  <- c / total
p_ApBp <- d / total

# Tabel peluang bersama
joint_prob <- matrix(c(p_AB, p_ABp, p_ApB, p_ApBp), nrow = 2, byrow = TRUE, dimnames = list(c("Ikut Pelatihan", "Tidak Ikut Pelatihan"), c("Lulus", "Tidak Lulus")))
joint_prob

##                          Lulus Tidak Lulus
## Ikut Pelatihan       0.4666667   0.1166667
## Tidak Ikut Pelatihan 0.1500000   0.2666667

Interpretasi:

P(Ikut Pelatihan \(\cap\) Lulus) = 0.467 → Artinya, 46.7% dari seluruh karyawan adalah peserta pelatihan yang lulus ujian.
P(Ikut Pelatihan \(\cap\) Tidak Lulus) = 0.117 → Sebanyak 11.7% dari seluruh karyawan ikut pelatihan namun tidak lulus.
P(Tidak Ikut Pelatihan \(\cap\) Lulus) = 0.15 → Hanya 15% dari seluruh karyawan yang tidak ikut pelatihan tapi lulus ujian.
P(Tidak Ikut Pelatihan \(\cap\) Tidak Lulus) = 0.267 → Terdapat sebanyak 26.7% karyawan yang tidak ikut pelatihan dan tidak lulus ujian.

- Peluang Marginal

Probabilitas dari satu variabel tanpa memperhatikan variabel lainnya, dihitung dengan menjumlahkan frekuensi baris atau kolom lalu dibagi total keseluruhan.

Rumus umum:

Peluang marginal baris:

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

Peluang marginal kolom:

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

Perhitungan Manual:

P(Ikut Pelatihan) \(= \frac{70}{120} = 0.583\)
P(Tidak Ikut Pelatihan) \(= \frac{50}{120} = 0.417\)
P(Lulus) \(= \frac{74}{120} = 0.617\)
P(Tidak Lulus) \(= \frac{46}{120} = 0.383\)

Perhitungan dengan R:

# Buat tabel kontingensi
tabel <- matrix(c(56, 14, 18, 32), nrow = 2, byrow = TRUE)
rownames(tabel) <- c("Ikut", "Tidak Ikut")
colnames(tabel) <- c("Lulus", "Tidak Lulus")
n <- sum(tabel)
tabel

##            Lulus Tidak Lulus
## Ikut          56          14
## Tidak Ikut    18          32

# Peluang marginal
p_A <- sum(tabel[1,]) / n
p_Ap <- sum(tabel[2,]) / n
p_B <- sum(tabel[,1]) / n
p_Bp <- sum(tabel[,2]) / n

marginal_prob <- list(
  "P(Ikut Pelatihan)" = p_A,
  "P(Tidak Ikut Pelatihan)" = p_Ap,
  "P(Lulus)" = p_B,
  "P(Tidak Lulus)" = p_Bp
)
marginal_prob

## $`P(Ikut Pelatihan)`
## [1] 0.5833333
## 
## $`P(Tidak Ikut Pelatihan)`
## [1] 0.4166667
## 
## $`P(Lulus)`
## [1] 0.6166667
## 
## $`P(Tidak Lulus)`
## [1] 0.3833333

Interpretasi:

P(Ikut Pelatihan) = 0.583 → Sebanyak 58.3% dari seluruh karyawan mengikuti pelatihan sebelum ujian.
P(Tidak Ikut Pelatihan) = 0.417 → Sebanyak 41.7% karyawan tidak mengikuti pelatihan.
P(Lulus) = 0.617 → Terdapat sebanyak 61.7% karyawan yang lulus ujian.
P(Tidak Lulus) = 0.383 → Sebanyak 38.3% karyawan tidak lulus ujian.

- Peluang Bersyarat (Conditional Probability)

Probabilitas suatu kejadian terjadi dengan syarat bahwa kejadian lain sudah terjadi.

Rumus umum:

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

Perhitungan Manual:

P(Lulus | Ikut Pelatihan) \(= \frac{56}{70} = 0.8\)
P(Lulus | Tidak Ikut Pelatihan) \(= \frac{18}{50} = 0.36\)

Perhitungan dengan R:

# Peluang bersyarat
p_B_given_A   <- p_AB / p_A
p_B_given_Ap  <- p_ApB / p_Ap

conditional_prob <- list(
  "P(Lulus | Ikut Pelatihan)" = p_B_given_A,
  "P(Lulus | Tidak Ikut Pelatihan)" = p_B_given_Ap
)
conditional_prob

## $`P(Lulus | Ikut Pelatihan)`
## [1] 0.8
## 
## $`P(Lulus | Tidak Ikut Pelatihan)`
## [1] 0.36

Interpretasi:

P(Lulus | Ikut Pelatihan) = 0.8 → Peluang karyawan lulus ujian jika ia mengikuti pelatihan adalah 80%. Ini menunjukkan bahwa pelatihan sangat berpengaruh terhadap kelulusan.
P(Lulus | Tidak Ikut Pelatihan) = 0.36 → Jika tidak ikut pelatihan, peluang lulus turun menjadi hanya 36%. Ini menguatkan dugaan bahwa pelatihan memberi dampak positif.

Ukuran Asosiasi

Ukuran asosiasi digunakan untuk mengukur kekuatan hubungan antara dua variabel kategorik dalam tabel kontingensi, terutama dalam konteks eksposur dan outcome.

Dalam kasus ini:

Eksposur: Mengikuti pelatihan (A)
Outcome: Lulus ujian (B)

- Beda Peluang (Risk Difference)

Beda peluang menunjukkan selisih peluang terjadinya suatu kejadian antara dua kelompok.

Rumus umum:

\[ RD = \left( \frac{n_{11}}{n_{1.}} \right) - \left( \frac{n_{21}}{n_{2.}} \right) = P(B|A) - P(B|A') \]

Perhitungan Manual:

\(RD = (\frac{56}{70}) - (\frac{18}{50}) = 0.44\)

Perhitungan dengan R:

risk_diff <- p_B_given_A - p_B_given_Ap
risk_diff

## [1] 0.44

Interpretasi:

RD = 0.44 → Peluang lulus meningkat sebesar 44% pada karyawan yang ikut pelatihan dibanding yang tidak ikut.

- Risiko Relatif (Relative Risk)

Risiko relatif membandingkan peluang terjadinya suatu kejadian antar dua kelompok. Nilai RR:

RR = 1 → tidak ada perbedaan risiko
RR > 1 → risiko lebih tinggi pada kelompok A
RR < 1 → risiko lebih rendah pada kelompok A

Rumus umum:

\[ RR = \frac{\frac{n_{11}}{n_{1.}}}{\frac{n_{21}}{n_{2.}}} = \frac{P(B \mid A)}{P(B \mid A')} \]

Perhitungan Manual:

\(RR = (\frac{56}{70})(\frac{50}{18}) = 2.22\)

Perhitungan dengan R:

relative_risk <- p_B_given_A / p_B_given_Ap
relative_risk

## [1] 2.222222

Interpretasi:

RR = 2.22 → Karyawan yang mengikuti pelatihan memiliki peluang 2.22 kali lebih besar untuk lulus ujian dibanding yang tidak ikut.

- Rasio Odds (Odds Ratio)

Odds Ratio (OR) membandingkan odds (peluang) terjadinya suatu kejadian antar dua kelompok. OR banyak digunakan dalam studi kasus-kontrol atau uji independensi.

Rumus umum:

\[ OR = \frac{n_{11} \times n_{22}}{n_{12} \times n_{21}} = \frac{P(B \mid A) / (1 - P(B \mid A))}{P(B \mid A') / (1 - P(B \mid A'))} \]

Perhitungan Manual:

\(OR = \frac{56\times32}{14\times18} = 7.11\)

Perhitungan dengan R:

odds_A  <- p_B_given_A / (1 - p_B_given_A)
odds_Ap <- p_B_given_Ap / (1 - p_B_given_Ap)
odds_ratio <- odds_A / odds_Ap
odds_ratio

## [1] 7.111111

Interpretasi:

OR = 7.11 → Odds (peluang relatif terhadap kegagalan) untuk lulus ujian bagi yang ikut pelatihan adalah 7.11 kali lebih tinggi dibanding yang tidak ikut pelatihan.

II. Inferensi Tabel Kontingensi Dua Arah

Inferensi dalam tabel kontingensi dua arah bertujuan untuk menentukan apakah terdapat hubungan yang signifikan secara statistik antara dua variabel kategorik. Artinya, kita ingin mengetahui apakah perbedaan proporsi yang tampak di tabel terjadi secara kebetulan atau memang mencerminkan hubungan nyata di populasi.

Untuk menguji hal ini, digunakan beberapa metode inferensial, seperti:

Uji Chi-Square (\(\chi^2\)) untuk independensi: Menguji apakah dua variabel saling bebas (tidak berhubungan).
Uji Fisher’s Exact Test: Alternatif saat ukuran sampel kecil atau ketika frekuensi harapan < 5.
Confidence interval untuk ukuran asosiasi (misalnya risk difference, relative risk, odds ratio): Memberikan kisaran nilai yang mungkin untuk parameter populasi.

Estimasi

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

- Estimasi Titik

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

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

di mana:

\(\hat{p}\) = estimasi titik proporsi
\(x\) = jumlah individu dalam kategori tertentu
\(n\) = total jumlah individu dalam sampel

- Estimasi Interval

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

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

di mana:

\(Z_{\alpha/2}\) = nilai dari distribusi normal standar untuk tingkat kepercayaan tertentu
\(\hat{p}\) = estimasi titik proporsi
\(n\) = ukuran sampel

Uji Hipotesis

Misalkan kita memiliki data berikut:

	Lulus	Tidak Lulus	Total
Ikut Pelatihan	56	14	70
Tidak Ikut Pelatihan	18	32	50
Total	74	46	120

- Uji Proporsi

Uji proporsi merupakan salah satu bentuk inferensi statistik yang digunakan untuk membandingkan proporsi kejadian antar dua kelompok.

Hipotesis

\(H_0\): \(p_1 = p_2\) (proporsi lulus sama pada kedua kelompok)
\(H_1\): \(p_1 \neq p_2\) (terdapat perbedaan proporsi lulus pada kedua kelompok)

Taraf Signifikansi

\(\alpha = 0.05\)

Statistik Uji

\[ z = \frac{p_1 - p_2}{\sqrt{p(1 - p) \left( \frac{1}{n_1} + \frac{1}{n_2} \right)}} \]

dengan:

\(p_1 = \frac{x_1}{n_1}\) = proporsi lulus kelompok ikut pelatihan
\(p_2 = \frac{x_2}{n_2}\) = proporsi lulus kelompok tidak ikut pelatihan
\(p = \frac{x_1 + x_2}{n_1 + n_2}\) = proporsi gabungan

Perhitungan:

# Jumlah data
x1 <- 56  # Lulus, ikut pelatihan
n1 <- 70  # Total ikut pelatihan
x2 <- 18  # Lulus, tidak ikut pelatihan
n2 <- 50  # Total tidak ikut pelatihan

# Proporsi masing-masing kelompok
p1 <- x1 / n1  # Proporsi kelulusan dari kelompok yang ikut pelatihan
p2 <- x2 / n2  # Proporsi kelulusan dari kelompok yang tidak ikut pelatihan

# Proporsi gabungan (pooled)
p <- (x1 + x2) / (n1 + n2)

# Statistik uji z
z_stat <- (p1 - p2) / sqrt(p * (1 - p) * (1/n1 + 1/n2))

# p-value (dua sisi)
p_value <- 2 * (1 - pnorm(abs(z_stat)))

# Output
list(
  "Proporsi lulus kelompok ikut pelatihan" = p1,
  "Proporsi lulus kelompok tidak ikut pelatihan" = p2,
  "Proporsi gabungan" = p,
  "z" = z_stat,
  "p-value" = p_value
)

## $`Proporsi lulus kelompok ikut pelatihan`
## [1] 0.8
## 
## $`Proporsi lulus kelompok tidak ikut pelatihan`
## [1] 0.36
## 
## $`Proporsi gabungan`
## [1] 0.6166667
## 
## $z
## [1] 4.887452
## 
## $`p-value`
## [1] 1.021492e-06

Kriteria Uji

Tolak \(H_0\) jika p-value < \(\alpha\), terima dalam hal lainnya.

Keputusan

Karena p-value (1.021492e-06) < \(\alpha\) (0.05), maka \(H_0\) ditolak.

Kesimpulan

Dengan tingkat kepercayaan sebesar 95%, dapat disimpulkan bahwa terdapat perbedaan proporsi kelulusan yang signifikan secara statistik antara karyawan yang mengikuti pelatihan dengan yang tidak. Dari perhitungan proporsi di atas, diperoleh hasil bahwa kelompok yang mengikuti pelatihan memiliki tingkat kelulusan yang jauh lebih tinggi (80%) dibanding kelompok yang tidak ikut pelatihan (36%). Dengan selisih proporsi sebesar 44%, mengikuti pelatihan berpengaruh kuat dan positif terhadap kelulusan peserta.

- Uji Asosiasi

Uji asosiasi dalam tabel kontingensi 2x2 bertujuan untuk mengukur hubungan antara dua variabel kategori. Tiga ukuran utama dalam uji asosiasi yaitu Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR).

Hipotesis

\(H_0\): Tidak ada asosiasi antara variabel pelatihan dan kelulusan
\(H_1\): Terdapat asosiasi antara variabel pelatihan dan kelulusan

Taraf Signifikansi

\(\alpha = 0.05\)

Statistik Uji

Statistik uji Z:

Risk Difference (RD)

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

Relative Risk (RR)

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

Odds Ratio (OR)

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

Perhitungan:

n11 <- 56; n12 <- 14; n21 <- 18; n22 <- 32
n1. <- n11 + n12; n2. <- n21 + n22

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

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

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

# Hasil
list(RD = rd, SE_RD = se_rd, Z_RD = z_rd, RR = rr, SE_Ln_RR = se_ln_rr, Z_RR = z_rr, OR = or, SE_Ln_OR = se_ln_or, Z_OR = z_or)

## $RD
## [1] 0.44
## 
## $SE_RD
## [1] 0.08302839
## 
## $Z_RD
## [1] 5.299392
## 
## $RR
## [1] 2.222222
## 
## $SE_Ln_RR
## [1] 0.1978054
## 
## $Z_RR
## [1] 4.036834
## 
## $OR
## [1] 7.111111
## 
## $SE_Ln_OR
## [1] 0.4196323
## 
## $Z_OR
## [1] 4.674708

Kriteria Uji

Tolak \(H_0\) jika \(|Z| > Z_{\alpha/2} = 1.96\), terima dalam hal lainnya.

Keputusan

Untuk RD: \(|Z_{RD}| = 5.30 > 1.96\) → Tolak \(H_0\)
Untuk RR: \(|Z_{RR}| = 4.04 > 1.96\) → Tolak \(H_0\)
Untuk OR: \(|Z_{OR}| = 4.67 > 1.96\) → Tolak \(H_0\)

Kesimpulan

Dengan tingkat kepercayaan sebesar 95%, dapat disimpulkan bahwa terdapat asosiasi yang signifikan secara statistik antara keikutsertaan dalam pelatihan dan kelulusan.

- Uji Independensi

Uji independensi adalah metode statistik yang digunakan untuk mengetahui apakah dua variabel kategorik saling bebas (tidak berhubungan) atau justru memiliki hubungan dalam suatu populasi.

Hipotesis

\(H_0\): Variabel pelatihan dan kelulusan bersifat independen (tidak ada hubungan)
\(H_1\): Variabel pelatihan dan kelulusan tidak bersifat independen (ada hubungan)

Taraf Signifikansi

\(\alpha = 0.05\)

Statistik Uji

Statistik chi-square dihitung dengan:

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

di mana:

\(O\): nilai observasi dalam tabel kontingensi
\(E\): nilai yang diharapkan
\(E_{ij} = \frac{(R_i \times C_j)}{N}\)

dengan:
- \(R_i\): total baris ke-\(i\)
- \(C_j\): total kolom ke-\(j\)
- \(N\): total sampel

Perhitungan:

# Data observasi
observed <- matrix(c(56, 14, 18, 32), 
                   nrow = 2, byrow = TRUE)

# Marginal total
row_totals <- rowSums(observed)
col_totals <- colSums(observed)
grand_total <- sum(observed)

# Hitung nilai harapan (expected)
expected <- outer(row_totals, col_totals) / grand_total

# Hitung nilai chi-square per sel
chi_square_cells <- (observed - expected)^2 / expected

# Total nilai chi-square
chi_square_stat <- sum(chi_square_cells)

# Derajat bebas (df) = (baris - 1)(kolom - 1) = 1
df <- (nrow(observed) - 1) * (ncol(observed) - 1)

# Hitung p-value
p_value <- pchisq(chi_square_stat, df = df, lower.tail = FALSE)

# Tampilkan hasil
list(
  "Chi-Square" = chi_square_stat,
  "p-value" = p_value
)

## $`Chi-Square`
## [1] 23.88719
## 
## $`p-value`
## [1] 1.021492e-06

Kriteria Uji

Tolak \(H_0\) jika p-value < \(\alpha\), terima dalam hal lainnya.

Keputusan

Karena p-value (1.021492e-06) < \(\alpha\) (0.05), maka \(H_0\) ditolak.

Kesimpulan

Dengan tingkat kepercayaan sebesar 95%, dapat disimpulkan bahwa terdapat hubungan yang signifikan secara statistik antara status pelatihan dan kelulusan ujian. Dengan kata lain, kemungkinan lulus ujian terkait dengan apakah seseorang mengikuti pelatihan atau tidak.

Analisis Residual dalam Tabel Kontingensi

- Analisis Pearson Residual

Misalkan kita memiliki data berikut:

	Lulus	Tidak Lulus	Total
Ikut Pelatihan	56	14	70
Tidak Ikut Pelatihan	18	32	50
Total	74	46	120

# Data Observasi
observed <- matrix(c(56, 14, 18, 32), nrow = 2, byrow = TRUE)

# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected

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

# Standardized Residual
row_sum <- rowSums(observed)
col_sum <- colSums(observed)
total_sum <- sum(observed)
standardized_residual <- (observed - expected) / 
  sqrt(expected * 
       outer(1 - row_sum / total_sum, 1 - col_sum / total_sum))

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

## $Pearson_Residual
##           [,1]      [,2]
## [1,]  1.953283 -2.477434
## [2,] -2.311156  2.931340
## 
## $Standardized_Residual
##           [,1]      [,2]
## [1,]  4.887452 -4.887452
## [2,] -4.887452  4.887452

Pearson Residual:

Nilai residu terbesar muncul di sel:
- (Ikut pelatihan, Tidak lulus) → −2.48
- (Tidak ikut pelatihan, Lulus) → −2.31
Artinya, jumlah aktual yang tidak lulus meskipun ikut pelatihan dan yang lulus tanpa ikut pelatihan lebih rendah dari yang diharapkan jika tidak ada hubungan.

Standardized Residual:

Semua nilai absolut standardized residual > 2 mengindikasikan penyimpangan signifikan dari ekspektasi.

Kesimpulan

Ada hubungan yang sangat signifikan antara ikut pelatihan dan kelulusan ujian. Pelatihan cenderung meningkatkan kemungkinan lulus.

- Deteksi Outlier Menggunakan Residual

# Data Observasi
observed <- matrix(c(56, 14, 18, 32), nrow = 2, byrow = TRUE)
colnames(observed) <- c("Lulus", "Tidak Lulus")
rownames(observed) <- c("Ikut Pelatihan", "Tidak Ikut Pelatihan")

# Hitung nilai ekspektasi
expected <- chisq.test(observed)$expected

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

# Standardized Residual
row_sum <- rowSums(observed)
col_sum <- colSums(observed)
total_sum <- sum(observed)
standardized_residual <- (observed - expected) / 
  sqrt(expected * 
       outer(1 - row_sum / total_sum, 1 - col_sum / total_sum))

# Menampilkan hasil
list(
Observed = observed,
Expected = expected,
Pearson_Residual = pearson_residual,
Standardized_Residual = standardized_residual
)

## $Observed
##                      Lulus Tidak Lulus
## Ikut Pelatihan          56          14
## Tidak Ikut Pelatihan    18          32
## 
## $Expected
##                         Lulus Tidak Lulus
## Ikut Pelatihan       43.16667    26.83333
## Tidak Ikut Pelatihan 30.83333    19.16667
## 
## $Pearson_Residual
##                          Lulus Tidak Lulus
## Ikut Pelatihan        1.953283   -2.477434
## Tidak Ikut Pelatihan -2.311156    2.931340
## 
## $Standardized_Residual
##                          Lulus Tidak Lulus
## Ikut Pelatihan        4.887452   -4.887452
## Tidak Ikut Pelatihan -4.887452    4.887452

Karena semua nilai residual |r| > 3, maka keempat sel merupakan outlier signifikan, menunjukkan hubungan yang sangat kuat antara keikutsertaan pelatihan dan kelulusan ujian.

III. Tabel Kontingensi Tiga Arah

Tabel kontingensi tiga arah adalah tabel yang digunakan untuk menyajikan data dari tiga variabel kategorik secara bersamaan. Tujuannya adalah untuk menganalisis hubungan antara dua variabel dengan mempertimbangkan pengaruh atau peran dari variabel ketiga, baik secara keseluruhan (marginal) maupun dalam tiap kelompok (parsial), serta menguji apakah hubungan tersebut konsisten di seluruh kelompok dengan bantuan uji seperti Cochran-Mantel-Haenszel.

Contoh Kasus

Akan diteliti hubungan antara kebiasaan merokok (Perokok/Non-Perokok) dan kejadian penyakit jantung (Ya/Tidak) dengan mempertimbangkan faktor usia (Muda/Tua) sebagai confounder. Data diperoleh dari hasil survei kesehatan yang dilakukan pada 400 responden yang disajikan dalam tabel kontingensi tiga arah sebagai berikut.

Usia	Merokok	Penyakit Jantung: Ya	Penyakit Jantung: Tidak	Total
Muda	Ya	30	70	100
Muda	Tidak	20	80	100
Tua	Ya	40	60	100
Tua	Tidak	25	75	100

Tabel Parsial & Tabel Marginal

- Tabel Parsial

Tabel parsial adalah tabel dua arah yang disusun untuk setiap kategori dari variabel ketiga, sehingga menunjukkan hubungan antara dua variabel utama dalam tiap kelompok.

Kelompok Usia Muda

Merokok	Penyakit Jantung: Ya	Penyakit Jantung: Tidak	Total
Ya	30	70	100
Tidak	20	80	100

Kelompok Usia Tua

Merokok	Penyakit Jantung: Ya	Penyakit Jantung: Tidak	Total
Ya	40	60	100
Tidak	25	75	100

- Tabel Marginal

Tabel marginal adalah tabel dua arah yang diperoleh dengan menggabungkan atau menjumlahkan data dari semua kategori variabel ketiga, sehingga menunjukkan hubungan dua variabel tanpa mempertimbangkan variabel ketiga.

Gabungan Usia Muda dan Tua

Merokok	Penyakit Jantung: Ya	Penyakit Jantung: Tidak	Total
Ya	70	130	200
Tidak	45	155	200

Distribusi Peluang

- Peluang Bersama (Joint Probability)

Peluang terjadinya tiga kejadian (dari tiga variabel kategorik) secara bersamaan.

Rumus umum:

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

dengan:

\(f(Z, X, Y)\) = jumlah observasi pada sel (Z, X, Y)
\(N\) = total jumlah observasi

Perhitungan:

data <- array(
  c(30, 20, 70, 80, 40, 25, 60, 75),
  dim = c(2, 2, 2),
  dimnames = list(
    Merokok = c("Ya", "Tidak"),
    Jantung = c("Ya", "Tidak"),
    Usia = c("Muda", "Tua")
  )
)
# Hitung peluang bersama
total_data <- sum(data)
joint_prob <- data / total_data
ftable(joint_prob)

##                 Usia   Muda    Tua
## Merokok Jantung                   
## Ya      Ya           0.0750 0.1000
##         Tidak        0.1750 0.1500
## Tidak   Ya           0.0500 0.0625
##         Tidak        0.2000 0.1875

Interpretasi:

Dari hasil perhitungan di atas, terlihat bahwa kebiasaan merokok dan kelompok usia berhubungan dengan peningkatan risiko penyakit jantung. Peluang terkena penyakit jantung lebih tinggi pada perokok dibandingkan non-perokok di setiap kelompok usia, serta lebih tinggi pada kelompok usia tua dibandingkan yang muda, baik perokok maupun non-perokok. Ini menunjukkan bahwa merokok dan bertambahnya usia berkontribusi terhadap peningkatan risiko penyakit jantung.

- Peluang Marginal

Peluang terjadinya satu atau dua variabel tanpa memperhatikan variabel lainnya, diperoleh dengan menjumlahkan peluang bersama pada variabel yang diabaikan.

Rumus umum:

\[ P(X) = \frac{\sum_{i} f_i}{n} \]

dengan:

\(P(X)\) = peluang marginal variabel X
\(f_i\) = frekuensi kategori tertentu
\(n\) = total observasi

Perhitungan:

marginal_usia <- apply(joint_prob, 3, sum)
marginal_merokok <- apply(joint_prob, 1, sum)
marginal_jantung <- apply(joint_prob, 2, sum)
list(
  "Peluang marginal usia" = marginal_usia,
  "Peluang marginal merokok" = marginal_merokok,
  "Peluang marginal jantung" = marginal_jantung
)

## $`Peluang marginal usia`
## Muda  Tua 
##  0.5  0.5 
## 
## $`Peluang marginal merokok`
##    Ya Tidak 
##   0.5   0.5 
## 
## $`Peluang marginal jantung`
##     Ya  Tidak 
## 0.2875 0.7125

Interpretasi:

Hasil perhitungan di atas menunjukkan bahwa proporsi individu dalam data terbagi rata antara kelompok usia muda dan tua (masing-masing 0.5), serta antara perokok dan non-perokok (masing-masing 0.5). Sementara itu, peluang marginal untuk penyakit jantung menunjukkan bahwa 28.75% dari total individu mengalami penyakit jantung, sedangkan 71.25% tidak, yang mengindikasikan bahwa mayoritas individu dalam data ini tidak memiliki penyakit jantung.

- Peluang Bersyarat (Conditional Probability)

Peluang suatu kejadian dari satu atau dua variabel dengan syarat nilai dari variabel lain diketahui.

Rumus umum:

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

dengan:

\(P(Z | X, Y)\) = peluang terjadinya Z jika X dan Y sudah terjadi
\(P(Z, X, Y)\) = peluang kejadian Z, X, dan Y terjadi bersama
\(P(X, Y)\) = peluang kejadian X dan Y

Perhitungan:

data <- matrix(c(30, 70, 20, 80, 40, 60, 25, 75), 
               nrow = 4, byrow = TRUE,
               dimnames = list(c("Muda, Merokok", "Muda, Tidak Merokok", "Tua, Merokok", "Tua, Tidak Merokok"), c("Jantung_Yes", "Jantung_No")))

total_data <- rowSums(data)  # Total per kategori

# Peluang bersyarat P(Jantung | Usia, Merokok)
prob_conditional <- data / total_data
prob_conditional

##                     Jantung_Yes Jantung_No
## Muda, Merokok              0.30       0.70
## Muda, Tidak Merokok        0.20       0.80
## Tua, Merokok               0.40       0.60
## Tua, Tidak Merokok         0.25       0.75

Interpretasi:

Hasil di atas menunjukkan bahwa perokok memiliki risiko lebih tinggi terkena penyakit jantung dibandingkan non-perokok di setiap kelompok usia. Pada kelompok usia muda, peluang terkena penyakit jantung adalah 30% bagi perokok dan 20% bagi non-perokok, sedangkan pada kelompok usia tua, peluangnya meningkat menjadi 40% bagi perokok dan 25% bagi non-perokok. Ini menunjukkan bahwa merokok meningkatkan risiko penyakit jantung, dan dampaknya semakin besar seiring bertambahnya usia.

Tabel Peluang Bersyarat

Peluang bersyarat dihitung berdasarkan jumlah total dalam setiap kategori

Untuk Usia Muda

Merokok	Jantung_Yes	Jantung_No
Merokok	0.30	0.70
Tidak Merokok	0.20	0.80

Untuk Usia Tua

Merokok	Jantung_Yes	Jantung_No
Merokok	0.40	0.60
Tidak Merokok	0.25	0.75

Ukuran Asosiasi

Ukuran asosiasi dalam tabel kontingensi tiga arah digunakan untuk mengukur kekuatan dan arah hubungan antara dua variabel dengan mempertimbangkan pengaruh variabel ketiga sebagai pengganggu. Dalam konteks ini, ukuran asosiasi dapat dihitung secara kondisional untuk setiap level variabel pengganggu, serta secara marginal untuk seluruh data. Perbandingan antara nilai-nilai ini membantu mengidentifikasi apakah terdapat confounding atau interaksi, serta apakah asosiasi antara dua variabel tetap konsisten di seluruh strata variabel ketiga.

- Beda Peluang (Risk Difference)

Beda peluang merupakan selisih peluang terjadinya suatu kejadian antara dua kelompok yang dihitung pada setiap level variabel ketiga untuk melihat pengaruh setelah dikendalikan.

Rumus umum:

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

dengan:

\(P(Z \mid X_1,Y)\) = probabilitas terjadinya penyakit jantung pada perokok di usia tertentu
\(P(Z \mid X_2,Y)\) = probabilitas terjadinya penyakit jantung pada non-perokok di usia tertentu

- Untuk usia muda

Perhitungan Manual:

\[ RD_{Muda} = (\frac{30}{100}) - (\frac{20}{100}) = 0.1 \]

Perhitungan dengan R:

RD_Muda <- (30/100) - (20/100)
RD_Muda

## [1] 0.1

Interpretasi:

Peluang terkena penyakit jantung pada perokok lebih tinggi sebesar 10% dibandingkan non-perokok pada kelompok usia muda.

- Untuk usia tua

Perhitungan Manual:

\[ RD_{Tua} = (\frac{40}{100}) - (\frac{25}{100}) = 0.15 \]

Perhitungan dengan R:

RD_Tua <- (40/100) - (25/100)
RD_Tua

## [1] 0.15

Interpretasi:

Peluang terkena penyakit jantung pada perokok lebih tinggi sebesar 15% dibandingkan non-perokok pada kelompok usia tua.

- Marginal

Perhitungan Manual:

\[ RD_{Marginal} = (\frac{70}{200}) - (\frac{45}{200}) = 0.125 \]

Perhitungan dengan R:

RD_Marginal <- (70/200) - (45/200)
RD_Marginal

## [1] 0.125

Interpretasi:

Secara keseluruhan, perokok memiliki peluang 12.5% lebih tinggi untuk terkena penyakit jantung dibandingkan non-perokok.

- Risiko Relatif (Relative Risk)

Risiko relatif membandingkan risiko kejadian antar dua kelompok pada setiap tingkat variabel ketiga. Ini membantu mengevaluasi apakah asosiasi antara dua variabel konsisten di seluruh strata variabel pengganggu.

Rumus umum:

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

dengan:

\(P(Z \mid X_1, Y)\) = probabilitas terjadinya penyakit jantung pada perokok di usia tertentu
\(P(Z \mid X_2, Y)\) = probabilitas terjadinya penyakit jantung pada non-perokok di usia tertentu

- Untuk usia muda

Perhitungan Manual:

\[ RR_{Muda} = (\frac{30}{100})(\frac{100}{20}) = 1.5 \]

Perhitungan dengan R:

RR_Muda <- (30/100) / (20/100)
RR_Muda

## [1] 1.5

Interpretasi:

Perokok memiliki risiko 1.5 kali lebih besar terkena penyakit jantung dibandingkan non-perokok pada kelompok usia muda.

- Untuk usia tua

Perhitungan Manual:

\[ RR_{Tua} = (\frac{40}{100})(\frac{100}{25}) = 1.6 \]

Perhitungan dengan R:

RR_Tua <- (40/100) / (25/100)
RR_Tua

## [1] 1.6

Interpretasi:

Perokok memiliki risiko 1.6 kali lebih besar terkena penyakit jantung dibandingkan non-perokok pada kelompok usia tua.

- Marginal:

Perhitungan Manual:

\[ RR_{Marginal} = (\frac{70}{200})(\frac{200}{45}) = 1.56 \]

Perhitungan dengan R:

RR_Marginal <- (70/200) / (45/200)
RR_Marginal

## [1] 1.555556

Interpretasi:

Secara keseluruhan, perokok memiliki risiko 1.56 kali lebih besar terkena penyakit jantung dibandingkan non-perokok.

- Rasio Odds (Odds Ratio)

Odds ratio mengukur kekuatan asosiasi antara dua variabel pada setiap tingkat variabel ketiga. Ini berguna untuk mengidentifikasi adanya interaksi atau efek perancu.

Rumus umum:

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

dengan:

\(P(Z \mid X_1, Y)\) = probabilitas terjadinya penyakit jantung pada perokok di usia tertentu
\(P(Z \mid X_2, Y)\) = probabilitas terjadinya penyakit jantung pada non-perokok di usia tertentu

- Untuk usia muda

Perhitungan Manual:

\[ OR_{Muda} = \frac{(30\times80)}{(20\times70)} = 1.71 \]

Perhitungan dengan R:

OR_Muda <- (30 * 80) / (20 * 70)
OR_Muda

## [1] 1.714286

Interpretasi:

Odds terkena penyakit jantung pada perokok 1.71 kali lebih besar dibandingkan non-perokok pada kelompok usia muda.

- Untuk usia tua

Perhitungan Manual:

\[ OR_{Tua} = \frac{(40\times75)}{(25\times60)} = 2 \]

Perhitungan dengan R:

OR_Tua <- (40 * 75) / (25 * 60)
OR_Tua

## [1] 2

Interpretasi:

Odds terkena penyakit jantung pada perokok 2 kali lebih besar dibandingkan non-perokok pada kelompok usia tua.

- Marginal

Perhitungan Manual:

\[ OR_{Marginal} = \frac{(70\times155)}{(45\times130)} = 1.85 \]

Perhitungan dengan R:

OR_Marginal <- (70 * 155) / (45 * 130)
OR_Marginal

## [1] 1.854701

Interpretasi:

Secara keseluruhan, odds terkena penyakit jantung pada perokok 1.85 kali lebih besar dibandingkan non-perokok.

Conditional Independence

Conditional independence (kemandirian bersyarat) dalam tabel kontingensi terjadi ketika dua variabel menjadi independen setelah dikendalikan oleh variabel ketiga.

Contoh Kasus

Misalkan kita memiliki data tentang tingkat stres (X), produktivitas kerja (Y), dan dukungan sosial (Z) dari karyawan suatu perusahaan.

X (Stres): Tinggi, Rendah
Y (Produktivitas): Tinggi, Rendah
Z (Dukungan Sosial): Ada, Tidak Ada

Tabel 1: Untuk Z = Ada (Dukungan Sosial Ada)

	Produktivitas Tinggi	Produktivitas Rendah
Stres Tinggi	20	10
Stres Rendah	40	30

Tabel 2: Untuk Z = Tidak Ada (Dukungan Sosial Tidak Ada)

	Produktivitas Tinggi	Produktivitas Rendah
Stres Tinggi	10	20
Stres Rendah	15	25

library(knitr)
library(kableExtra)

# Data array 3 dimensi: (X = Stres, Y = Produktivitas, Z = Dukungan Sosial)
stress_productivity <- array(
  c(20, 40, 10, 30,  # Z = Ada (dimensi ke-3 = 1)
    10, 15, 20, 25), # Z = Tidak Ada (dimensi ke-3 = 2)
  dim = c(2, 2, 2),
  dimnames = list(
    "Stres" = c("Tinggi", "Rendah"),
    "Produktivitas" = c("Tinggi", "Rendah"),
    "Dukungan" = c("Ada", "Tidak Ada")
  )
)

# Menampilkan tabel per kondisi Dukungan Sosial
for (z in dimnames(stress_productivity)$Dukungan) {
  df <- as.data.frame.matrix(stress_productivity[,,z])
  colnames(df) <- c("Tinggi", "Rendah")
  rownames(df) <- c("Stres Tinggi", "Stres Rendah")
  print(kable(df, format = "html", booktabs = TRUE,
              caption = paste("Tabel untuk Dukungan Sosial =", z)) %>%
          kable_styling(full_width = FALSE, position = "center"))
}

## <table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
## <caption>Tabel untuk Dukungan Sosial = Ada</caption>
##  <thead>
##   <tr>
##    <th style="text-align:left;">   </th>
##    <th style="text-align:right;"> Tinggi </th>
##    <th style="text-align:right;"> Rendah </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;"> Stres Tinggi </td>
##    <td style="text-align:right;"> 20 </td>
##    <td style="text-align:right;"> 10 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Stres Rendah </td>
##    <td style="text-align:right;"> 40 </td>
##    <td style="text-align:right;"> 30 </td>
##   </tr>
## </tbody>
## </table><table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
## <caption>Tabel untuk Dukungan Sosial = Tidak Ada</caption>
##  <thead>
##   <tr>
##    <th style="text-align:left;">   </th>
##    <th style="text-align:right;"> Tinggi </th>
##    <th style="text-align:right;"> Rendah </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;"> Stres Tinggi </td>
##    <td style="text-align:right;"> 10 </td>
##    <td style="text-align:right;"> 20 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Stres Rendah </td>
##    <td style="text-align:right;"> 15 </td>
##    <td style="text-align:right;"> 25 </td>
##   </tr>
## </tbody>
## </table>

Pengujian Conditional Independence

# Uji Chi-Square untuk masing-masing dukungan sosial
chisq_ada_dukungan <- chisq.test(stress_productivity[,,"Ada"])
chisq_tidak_dukungan <- chisq.test(stress_productivity[,,"Tidak Ada"])

chisq_ada_dukungan

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  stress_productivity[, , "Ada"]
## X-squared = 0.44643, df = 1, p-value = 0.504

chisq_tidak_dukungan

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  stress_productivity[, , "Tidak Ada"]
## X-squared = 0.011667, df = 1, p-value = 0.914

Interpretasi Hasil

Karena p-value keduanya > 0.05, maka tidak ada hubungan signifikan antara stres dan produktivitas setelah dikendalikan oleh ada dan tidaknya dukungan sosial.

IV. Inferensi Tabel Kontingensi Tiga Arah

Inferensi pada tabel kontingensi tiga arah digunakan untuk memahami hubungan antara dua variabel utama sambil mengendalikan variabel ketiga sebagai perancu atau efek interaksi. Dengan analisis ini, kita bisa menguji apakah asosiasi bersifat homogen di seluruh strata atau menghitung ukuran gabungan saat tidak ada interaksi signifikan.

Pengujian Independensi

Pengujian independensi dalam tabel kontingensi tiga arah adalah metode statistik untuk menguji apakah dua variabel saling bebas (tidak berasosiasi) setelah memperhitungkan pengaruh variabel ketiga. Tujuan utamanya adalah mengevaluasi apakah hubungan antara dua variabel utama tetap ada atau hilang ketika stratifikasi dilakukan berdasarkan variabel tambahan tersebut. Salah satu metode yang umum digunakan adalah uji Cochran-Mantel-Haenszel (CMH), yang secara khusus dirancang untuk menguji asosiasi dua variabel dengan mengontrol satu variabel pengganggu.

Metode Cochran-Mantel-Haenszel (CMH)

Metode Cochran-Mantel-Haenszel (CMH) digunakan untuk menguji independensi antara dua variabel setelah mengendalikan variabel ketiga (strata). CMH menggabungkan informasi dari semua strata untuk menguji apakah ada asosiasi yang konsisten antara dua variabel utama, dengan asumsi bahwa odds ratio (atau relative risk) homogen di seluruh strata.

Misalkan kita memiliki data sebagai berikut:

Usia	Merokok	Penyakit Jantung: Ya	Penyakit Jantung: Tidak	Total
Muda	Ya	30	70	100
Muda	Tidak	20	80	100
Tua	Ya	40	60	100
Tua	Tidak	25	75	100

Hipotesis

\(H_0: \theta_{xy(k)} = 1 \quad \text{untuk setiap } k = 1, 2, \dots, K\) (Tidak ada asosiasi antara merokok dan penyakit jantung setelah mengontrol faktor usia)
\(H_1: \theta_{xy(k)} \ne 1 \quad \text{untuk paling sedikit satu } k\) (Ada asosiasi antara merokok dan penyakit jantung setelah mengontrol faktor usia)

Taraf Signifikansi

\(\alpha = 0.05\)

Statistik Uji

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

dengan:

\(n_{11k}\): nilai frekuensi sel baris 1 kolom 1 pada tabel parsial ke-\(k\)
\(\mu_{11k}\): nilai ekspektasi sel baris 1 kolom 1 pada tabel parsial ke-\(k\), dihitung dengan rumus:

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

Varians dari \(n_{11k}\) dihitung dengan rumus:

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

CMH mengikuti distribusi Chi-square dengan db = 1

Perhitungan:

datarray <- array(
  c(30, 20, 70, 80, 40, 25, 60, 75),
  dim = c(2, 2, 2),
  dimnames = list(
    Merokok = c("Ya", "Tidak"),
    Jantung = c("Ya", "Tidak"),
    Usia = c("Muda", "Tua")
  )
)
cmh_base <- mantelhaen.test(datarray, correct = FALSE)
cmh_base

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  datarray
## Mantel-Haenszel X-squared = 7.6421, df = 1, p-value = 0.005702
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.196221 2.898546
## sample estimates:
## common odds ratio 
##          1.862069

Nilai Kritis

chi_critical <- qchisq(0.95, df = 1)
chi_critical

## [1] 3.841459

Kriteria Uji

Tolak \(H_0\) jika \(CMH > \chi^2_{(1)}\) dan p-value < \(\alpha\), terima dalam hal lainnya.

Keputusan

Karena \(CMH\) (7.6421) > \(\chi^2_{(1)}\) (3.841459) dan p-value (0.005702) < \(\alpha\) (0.05), maka \(H_0\) ditolak.

Kesimpulan

Hasil uji Cochran-Mantel-Haenszel menunjukkan bahwa setelah mengontrol faktor usia, terdapat hubungan yang signifikan antara merokok dan penyakit jantung. Hal ini menunjukkan bahwa merokok merupakan faktor risiko independen yang berkontribusi terhadap penyakit jantung.

V. Generalized Linear Model (GLM)

Generalized Linear Model (GLM) adalah perluasan dari model regresi linear yang memungkinkan kita memodelkan hubungan antara satu atau lebih variabel prediktor dengan variabel respons yang tidak harus berdistribusi normal. GLM terdiri dari tiga komponen utama:

Distribusi: Variabel respons \(Y\) berasal dari keluarga distribusi eksponensial umum (contoh: normal, binomial, poisson).
Fungsi link: Menghubungkan nilai harapan dari \(Y\), yaitu \(\mu = E(Y)\), dengan prediktor linear melalui fungsi link:

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

Fungsi sistematik (linear predictor): Kombinasi linear dari variabel prediktor:

\[ \eta = \beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p \]

Exponential Family (Keluarga Eksponensial)

Keluarga eksponensial adalah kelas umum dari distribusi probabilitas yang dapat dituliskan dalam bentuk tertentu yang melibatkan eksponensial dan memiliki bentuk umum:

\[ f(y; \theta) = s(y) \, t(\theta) \, \exp \left( a(y) b(\theta) \right) \]

di mana \(a, b, s\) dan \(t\) merupakan fungsi yang diketahui dan \(\theta\) merupakan parameter.

Contoh

Distribusi Binomial

Untuk \(Y \sim \text{Binomial}(n, p)\), bentuk eksponensialnya:

\[ f_Y(y) = \binom{n}{y} p^y (1-p)^{n-y} = \binom{n}{y} \exp \left[ y \log \left( \frac{p}{1 - p} \right) + n \log(1 - p) \right] \]

Distribusi Poisson

\[ f_Y(y; \lambda) = \frac{e^{-\lambda} \lambda^y}{y!} = \frac{1}{y!} \exp \left( y \log \lambda - \lambda \right) \]

Distribusi Normal (dengan varians diketahui)

\[ f_Y(y; \mu) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( -\frac{(y - \mu)^2}{2\sigma^2} \right) \]

Model Regresi Logistik

Model regresi logistik adalah salah satu jenis model dalam Generalized Linear Model (GLM), yang digunakan untuk memodelkan variabel dependen kategorikal dengan dua kategori (disebut sebagai biner) atau lebih. Model ini sering digunakan ketika variabel dependen merupakan hasil klasifikasi, misalnya “ya” atau “tidak”, “sukses” atau “gagal”.

Fungsi link yang digunakan dalam regresi logistik adalah logit, yang merupakan logaritma dari rasio peluang (odds ratio).

Fungsi logit adalah:

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

Regresi logistik biasanya digunakan untuk memprediksi probabilitas suatu kejadian berdasarkan satu atau lebih variabel prediktor. Model regresi logistik dapat ditulis sebagai:

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

di mana:

Setelah menghitung nilai logit, kita dapat menghitung probabilitas kejadian dengan mengubah logit kembali menggunakan fungsi sigmoid atau logistik:

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

Fungsi sigmoid ini menghasilkan nilai antara 0 dan 1, yang sesuai untuk merepresentasikan probabilitas.

Contoh

Simulasi Data

Kita akan mensimulasikan data untuk sebuah studi pengaruh kadar kolesterol terhadap risiko penyakit jantung. - Variabel prediktor: kadar kolesterol (data kontinu) - Variabel dependen: penyakit jantung (data biner: 1 = ya, 0 = tidak)

set.seed(42)
n <- 300
X <- rnorm(n, mean = 200, sd = 30)
beta_0 <- -8
beta_1 <- 0.03
log_odds <- beta_0 + beta_1 * X
p <- 1 / (1 + exp(-log_odds))
Y <- rbinom(n, 1, p)
sim.data <- data.frame(X, Y)
head(sim.data)

Estimasi Model Regresi Logistik

model <- glm(Y ~ X, data = sim.data, family = "binomial")
summary(model)

## 
## Call:
## glm(formula = Y ~ X, family = "binomial", data = sim.data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -8.274495   1.366315  -6.056 1.39e-09 ***
## X            0.032084   0.006321   5.076 3.85e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 270.34  on 299  degrees of freedom
## Residual deviance: 239.46  on 298  degrees of freedom
## AIC: 243.46
## 
## Number of Fisher Scoring iterations: 5

Interpretasi Koefisien

Untuk menginterpretasikan koefisien, hitung odds ratio (OR) yang merupakan eksponen dari koefisien.

exp(coef(model))

##  (Intercept)            X 
## 0.0002549367 1.0326046464

Dari hasil di atas, diperoleh odds ratio untuk X adalah 1.033. Artinya, setiap kenaikan satu unit pada X meningkatkan peluang kejadian Y = 1 sebanyak 1.033 kali.

Prediksi dan Visualisasi Kurva Logit

Memprediksi probabilitas kejadian Y = 1 berdasarkan nilai X dan menggambarkan kurva logit.

# Prediksi probabilitas
sim.data$predicted_prob <- predict(model, type = "response")

# Visualisasi kurva logit
library(ggplot2)
ggplot(sim.data, aes(x = X, y = predicted_prob)) +
  geom_point(aes(color = factor(Y)), alpha = 0.5) +
  geom_line(aes(y = predicted_prob), color = "red") +
  labs(title = "Kurva Logit dan Prediksi Probabilitas", x = "Variabel Prediktor X", y = "Prediksi Probabilitas P(Y=1)") +
  theme_minimal()

Interpretasi:

Hasil grafik menunjukkan bahwa semakin tinggi kadar kolesterol (X), semakin besar kemungkinan seseorang mengidap penyakit jantung (Y = 1). Model logistik berhasil membedakan dua kelas (Y = 0 dan Y = 1) dengan baik, mengikuti pola sigmoid yang wajar.

Evaluasi Model

# Membuat prediksi biner (0 atau 1)
sim.data$predicted_class <- ifelse(sim.data$predicted_prob > 0.5, 1, 0)

# Menghitung akurasi
accuracy <- mean(sim.data$predicted_class == sim.data$Y)
accuracy

## [1] 0.8266667

# Confusion Matrix
table(Predicted = sim.data$predicted_class, Actual = sim.data$Y)

##          Actual
## Predicted   0   1
##         0 245  47
##         1   5   3

Interpretasi:

Model regresi logistik yang dibuat berhasil mengklasifikasikan data dengan akurasi sebesar 82.7%, artinya sekitar 83 dari 100 orang diklasifikasikan dengan benar apakah mengidap penyakit jantung atau tidak.

Berdasarkan confusion matrix:

Model sangat baik dalam mengenali individu yang tidak mengidap penyakit jantung, dengan 245 prediksi benar dari 250 kasus aktual.
Namun, model kurang mampu mengenali individu yang mengidap penyakit jantung, karena hanya 3 dari 50 kasus aktual yang berhasil diprediksi dengan benar.

Model Regresi Poisson

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

Distribusi Poisson memiliki fungsi probabilitas:

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

Bentuk dalam format exponential family:

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

dengan:

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

Contoh

Misalkan terdapat dataset simulasi yang merepresentasikan jumlah pelanggan yang datang ke sebuah toko setiap hari berdasarkan dua variabel prediktor: jumlah iklan dan hari kerja (weekday atau bukan).

Simulasi Data

set.seed(123)
n <- 100
data.sim <- data.frame(
  ads = rpois(n, lambda = 5),
  weekday = sample(c(0, 1), size = n, replace = TRUE),
  customers = NA
)
lambda <- exp(1 + 0.2 * data.sim$ads + 0.5 * data.sim$weekday)
data.sim$customers <- rpois(n, lambda = lambda)
head(data.sim)

Variabel:

ads: jumlah iklan yang ditayangkan
weekday: apakah hari tersebut hari kerja (1) atau akhir pekan (0)
customers: jumlah pelanggan yang datang (data cacah)

Estimasi Model Regresi Poisson

Akan dilakukan pemodelan jumlah pelanggan (customers) menggunakan regresi Poisson dengan dua prediktor: ads dan weekday.

model <- glm(customers ~ ads + weekday, data = data.sim, family = poisson(link = "log"))
summary(model)

## 
## Call:
## glm(formula = customers ~ ads + weekday, family = poisson(link = "log"), 
##     data = data.sim)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.99135    0.09825  10.090   <2e-16 ***
## ads          0.19490    0.01360  14.330   <2e-16 ***
## weekday      0.53054    0.06425   8.258   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 359.28  on 99  degrees of freedom
## Residual deviance: 112.45  on 97  degrees of freedom
## AIC: 523.29
## 
## Number of Fisher Scoring iterations: 4

Interpretasi Koefisien

exp(coef(model))

## (Intercept)         ads     weekday 
##    2.694866    1.215191    1.699852

Intercept = 2.69

Saat tidak ada iklan dan bukan merupakan hari kerja, diperkirakan ada sekitar 3 pelanggan yang datang (dibulatkan).
ads = 1.215

Setiap tambahan 1 unit iklan diperkirakan meningkatkan rata-rata jumlah pelanggan sebesar 21.5%, jika faktor lainnya tetap.
weekday = 1.7

Pada hari kerja (weekday = 1), rata-rata jumlah pelanggan 1.7 kali lebih banyak dibandingkan akhir pekan, dengan asumsi jumlah iklan tetap.

Visualisasi Prediksi

Memprediksi dan memvisualisasikan rata-rata pelanggan terhadap jumlah iklan, untuk hari kerja dan akhir pekan.

library(ggplot2)

# Buat grid prediksi
newdata <- expand.grid(
  ads = 0:10,
  weekday = c(0, 1)
)

# Tambahkan prediksi rata-rata (lambda)
newdata$predicted_customers <- predict(model, newdata = newdata, type = "response")

# Plot
ggplot(newdata, aes(x = ads, y = predicted_customers, color = factor(weekday))) +
  geom_line(size = 1.2) +
  labs(x = "Jumlah Iklan", y = "Prediksi Pelanggan",
       color = "Hari Kerja",
       title = "Prediksi Pelanggan Berdasarkan Jumlah Iklan dan Hari") +
  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.

Interpretasi:

Dari hasil grafik tersebut diperoleh informasi bahwa semakin banyak iklan, jumlah pelanggan akan semakin meningkat dan peningkatannya bersifat eksponensial. Selain itu, garis berwarna biru (weekday) selalu lebih tinggi dibandingkan garis berwarna merah (weekend), yang berarti jumlah pelanggan lebih banyak saat hari kerja dibandingkan akhir pekan untuk jumlah iklan yang sama. Jadi, dapat disimpulkan bahwa jumlah iklan dan hari kerja sama-sama mempunyai pengaruh positif terhadap jumlah pelanggan, dan efeknya saling memperkuat.

Evaluasi Model

Untuk mengevaluasi kesesuaian model, dapat dilakukan dengan melihat pola residual. Residual yang baik seharusnya tersebar secara acak di sekitar nol tanpa pola tertentu.

# Buat data frame residual
res_data <- data.frame(
  Index = 1:length(model$residuals),
  Residual = model$residuals
)

# Plot residual
ggplot(res_data, aes(x = Index, y = Residual)) +
  geom_point(color = "blue") +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residual Plot", x = "Index", y = "Residual")

Interpretasi:

Hasil plot residual di atas menunjukkan bahwa titik-titik residual tersebar secara acak di sekitar garis horizontal nol tanpa membentuk suatu pola tertentu. Ini mengindikasikan bahwa model sudah cukup baik.

VI. Inferensi Generalized Linear Model (GLM)

Inferensi dalam Generalized Linear Model (GLM) adalah proses menarik kesimpulan tentang hubungan antara variabel prediktor dan respon, dengan mempertimbangkan bahwa respon bisa berasal dari distribusi selain normal (seperti binomial atau Poisson). Tujuannya untuk menguji pengaruh variabel, memperkirakan parameter, dan menilai kecocokan model.

Ekspektasi dan Varians dalam GLM

Ekspektasi Estimator

Ekspektasi menunjukkan apakah suatu estimator tak bias, yaitu:

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

Dalam GLM, MLE bersifat asymptotically unbiased.

Varians Estimator

Varians menunjukkan presisi dari estimasi parameter:

\[ \mathrm{Var}(\hat{\beta}) \approx \left[ \mathbf{X}^\top \mathbf{W} \mathbf{X} \right]^{-1} \]

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

Distribusi Asimptotik Estimator

Dengan ukuran sampel besar:

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

Distribusi ini adalah dasar dari:

Uji Wald
Confdence interval
P-value

Varians dalam GLM Tidak Konstan

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

dengan:

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

Contoh:

Poisson: \(V(\mu) = \mu\)
Binomial: \(V(\mu) = \mu(1 - \mu)\)

Inferensi Parameter

Contoh Kasus

Sebuah penelitian ingin mengetahui apakah lama belajar (dalam jam) memengaruhi kemungkinan seorang siswa lulus ujian sertifikasi.

Data dikumpulkan dari 100 siswa dengan variabel:

x: lama belajar per minggu
y: status kelulusan (1 = lulus, 0 = tidak lulus)

Model yang digunakan adalah regresi logistik.

- Uji Wald

Uji Wald digunakan untuk menguji signifikansi koefisien regresi. Tujuannya adalah mengetahui apakah setiap parameter \(\beta_j\) secara statistik berbeda dari nol (atau nilai tertentu).

Hipotesis

\(H_0\): \(\beta_j = 0\) (lama belajar tidak berpengaruh terhadap peluang lulus siswa)
\(H_1\): \(\beta_j \neq 0\) (lama belajar berpengaruh terhadap peluang lulus siswa)

Taraf Signifikansi

\(\alpha = 0.05\)

Statistik Uji

\[ W = Z^2 = (\frac{\hat\beta_j}{SE(\hat\beta_j)})^2 \]

Perhitungan:

set.seed(123)
n <- 100
x <- rnorm(n, mean = 5, sd = 2) # rata-rata belajar 5 jam/minggu
log_odds <- -2 + 0.8 * x
p <- 1 / (1 + exp(-log_odds))
y <- rbinom(n, 1, p)
data <- data.frame(x, y)

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

## 
## Call:
## glm(formula = y ~ x, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -4.7341     1.2962  -3.652  0.00026 ***
## x             1.3411     0.3088   4.343  1.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: 107.855  on 99  degrees of freedom
## Residual deviance:  66.221  on 98  degrees of freedom
## AIC: 70.221
## 
## Number of Fisher Scoring iterations: 6

Ambil nilai koefisien dan SE

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

Hitung statistik Z

Z <- beta_hat / se_beta
Z

##        x 
## 4.343132

Hitung statistik Wald

Wald_stat <- Z^2
Wald_stat

##       x 
## 18.8628

Hitung p-value

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

##            x 
## 1.404654e-05

Kriteria Uji

Tolak \(H_0\) jika p-value < \(\alpha\), terima dalam hal lainnya.

Keputusan

Karena p-value (1.404654e-05) < \(\alpha\) (0.05), maka \(H_0\) ditolak.

Kesimpulan

Dengan tingkat kepercayaan sebesar 95%, dapat disimpulkan bahwa lama belajar berpengaruh signifikan terhadap status kelulusan siswa.

- Uji Likelihood Ratio (Chi-Square)

Uji ini mengukur apakah penambahan variabel x secara signifikan meningkatkan kecocokan model terhadap data atau tidak.

Hipotesis

\(H_0\): Model tanpa prediktor sudah cukup baik (koefisien x = 0)
\(H_1\): Model dengan prediktor secara signifikan lebih baik

Pengujian

# Model 1: hanya intercept (model sederhana / null model)
model_null <- glm(y ~ 1, data = data, family = binomial)

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

Evaluasi Kebaikan Model

Akaike Information Criterion (AIC)

AIC(model_null, model)

Bayesian Information Criterion (BIC)

BIC(model_null, model)

Kesimpulan

Dari hasil uji di atas, diperoleh p-value (1.1e-10) < 0.05, yang berarti tolak \(H_0\) (model dengan prediktor secara signifikan lebih baik). Selain itu, berdasarkan hasil evaluasi kebaikan model menggunakan AIC dan BIC, diperoleh nilai AIC dan BIC yang lebih kecil pada model penuh dibandingkan model null sehingga dapat disimpulkan bahwa penambahan variabel prediktor (lama belajar) meningkatkan kualitas model secara signifikan.

VII. Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio

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

Nominal: Variabel yang tidak memiliki urutan logis antar kategorinya, hanya sebagai pembeda.

Contoh: jenis kelamin (laki-laki, perempuan). Dalam analisis regresi logistik, data nominal diubah menjadi variabel dummy (misal: Female = 1, Male = 0)

Ordinal: Variabel yang memiliki urutan logis antar kategori, tetapi jarak antar kategori belum tentu sama.

Contoh: tingkat pendidikan (SMA < Sarjana < Master < Doktor). Dalam analisis, variabel ordinal bisa diperlakukan:

– Sebagai nominal dengan membuat dummy untuk setiap kategori.

– Atau diperlakukan seperti rasio dengan memberi nilai tertentu (misalnya SMA = 1, Sarjana = 2, Master = 3, Doktor = 4) untuk mencerminkan urutan dan jarak antar tingkatannya.

Rasio: Variabel numerik kontinu yang memiliki nol absolut dan rasio bermakna. Contoh: pendapatan (dalam juta rupiah per bulan). Data ini dapat langsung digunakan dalam model tanpa transformasi khusus.

Contoh Kasus

Misalkan terdapat dataset simulasi yang merepresentasikan 500 karyawan pada sebuah perusahaan. Variabel yang diamati meliputi:

department (nominal): departemen tempat bekerja (HR, IT, Sales)
job_level (ordinal): level jabatan (Junior, Mid, Senior)
work_hours (rasio): jumlah jam kerja per minggu
promoted (biner): apakah karyawan tersebut mendapatkan promosi (1 = ya, 0 = tidak)

Kita ingin mengetahui pengaruh departemen, level jabatan, dan jam kerja terhadap kemungkinan karyawan mendapatkan promosi menggunakan regresi logistik.

Simulasi Data

n <- 500
set.seed(123)

department <- sample(c("HR", "IT", "Sales"), n, replace = TRUE)
job_level <- sample(c("Junior", "Mid", "Senior"), n, replace = TRUE, prob = c(0.5, 0.3, 0.2))
work_hours <- rnorm(n, mean = 40, sd = 5)

logit_p <- -3 + 
  0.6 * (department == "IT") + 
  1.0 * (department == "Sales") + 
  0.9 * as.numeric(factor(job_level, levels = c("Junior", "Mid", "Senior"), ordered = TRUE)) +
  0.07 * work_hours

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

promoted <- rbinom(n, 1, p)

sim_data <- tibble(promoted, department, job_level, work_hours)

head(sim_data)

Eksplorasi Data

sim_data %>%
  dplyr::group_by(promoted) %>%
  dplyr::summarise(
    Jumlah = dplyr::n(),
    Rata2_Work_Hours = mean(work_hours)
)

Berdasarkan hasil di atas, terlihat bahwa dari total 500 karyawan, 92 orang tidak dipromosikan dan 408 orang dipromosikan. Rata-rata jam kerja per minggu bagi yang tidak dipromosikan adalah sekitar 39.38 jam, sedangkan bagi yang dipromosikan adalah sekitar 40.25 jam. Artinya, karyawan yang mendapatkan promosi cenderung memiliki jam kerja mingguan sedikit lebih tinggi dibandingkan yang tidak dipromosikan. Ini mengindikasikan bahwa jam kerja mungkin berpengaruh terhadap peluang promosi.

Perlakuan Variabel Ordinal

Treat sebagai Nominal (Dummy)

library(dplyr)

## 
## Attaching package: 'dplyr'

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

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

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

sim_data_nominal <- sim_data %>%
  mutate(
    job_level = factor(job_level, levels = c("Junior", "Mid", "Senior"))
  )

model_nominal <- glm(promoted ~ department + job_level + work_hours, data = sim_data_nominal, family = binomial)
summary(model_nominal)

## 
## Call:
## glm(formula = promoted ~ department + job_level + work_hours, 
##     family = binomial, data = sim_data_nominal)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -0.83839    0.95669  -0.876  0.38085    
## departmentIT     0.75280    0.27823   2.706  0.00682 ** 
## departmentSales  1.55754    0.34914   4.461 8.15e-06 ***
## job_levelMid     1.20078    0.29219   4.110 3.96e-05 ***
## job_levelSenior  2.85784    0.73402   3.893 9.89e-05 ***
## work_hours       0.02724    0.02360   1.155  0.24827    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 477.40  on 499  degrees of freedom
## Residual deviance: 407.53  on 494  degrees of freedom
## AIC: 419.53
## 
## Number of Fisher Scoring iterations: 6

Interpretasi Koefisien

Intercept (-0.83839)
- Ini adalah log-odds dasar untuk individu pada level pekerjaan Junior (baseline), di departemen lain selain IT dan Sales, dan dengan work_hours = 0.
- Tidak signifikan (p = 0.38085), artinya baseline ini tidak berbeda secara signifikan dari probabilitas 50%.
- Odds Ratio = 0.43: Peluang dipromosikan lebih kecil daripada baseline (Junior level), namun tidak signifikan.
- Saat tidak ada iklan dan bukan merupakan hari kerja, diperkirakan ada sekitar 3 pelanggan yang datang (dibulatkan).
departmentIT (+0.75280)
- Individu di departemen IT memiliki log-odds promosi lebih tinggi sebesar 0.75 dibandingkan dengan departemen referensi (non-IT dan non-Sales).
- p = 0.00682 (signifikan di 1%), menunjukkan bahwa bekerja di departemen IT secara signifikan meningkatkan peluang promosi.
- Odds Ratio = 2.12: Peluang promosi di departemen IT sekitar 2 kali lebih besar dibandingkan referensi.
departmentSales (+1.55754)
- Individu di departemen Sales memiliki log-odds promosi lebih tinggi sebesar 1.56 dibandingkan departemen referensi.
- p < 0.0001 (sangat signifikan), artinya departemen Sales sangat berpengaruh terhadap peluang promosi.
- Odds Ratio = 4.74: Peluang promosi di departemen Sales sekitar 4.7 kali lebih besar dari referensi.
job_levelMid (+1.20078)
- Individu dengan level pekerjaan Mid memiliki log-odds promosi lebih tinggi sebesar 1.20 dibandingkan dengan Junior.
- p < 0.0001 (sangat signifikan), menunjukkan bahwa level pekerjaan Mid signifikan meningkatkan peluang promosi dibandingkan Junior.
- Odds Ratio = 3.32: Peluang promosi sekitar 3.3 kali lebih besar dari Junior.
job_levelSenior (+2.85784)
- Individu di level pekerjaan Senior memiliki log-odds promosi lebih tinggi sebesar 2.86 dibandingkan Junior.
- p < 0.0001 (sangat signifikan), artinya Senior level secara signifikan meningkatkan peluang promosi.
- Odds Ratio = 17.43: Peluang promosi pada level Senior sekitar 17 kali lebih besar dibandingkan Junior.
work_hours (+0.02724)
- Setiap penambahan 1 jam kerja dikaitkan dengan peningkatan log-odds promosi sebesar 0.027.
- Tidak signifikan (p = 0.24827), artinya jumlah jam kerja belum terbukti secara statistik memengaruhi peluang promosi.
- Odds Ratio = 1.03: Setiap tambahan 1 jam kerja meningkatkan peluang promosi sekitar 3%, namun tidak signifikan.

Interpretasi Goodness-of-Fit

Null deviance (477.40): Deviasi model tanpa prediktor.
Residual deviance (407.53): Deviasi model setelah memasukkan prediktor.
Penurunan nilai deviance menunjukkan bahwa model dengan prediktor membawa informasi yang cukup baik untuk menjelaskan variabel respons.
AIC (419.53): Ukuran untuk membandingkan model; semakin kecil AIC, semakin baik keseimbangan antara goodness-of-fit dan kompleksitas model.

Signifikansi Model

Variabel department, terutama Sales dan IT, secara signifikan meningkatkan peluang promosi.
Level pekerjaan Mid dan Senior juga sangat signifikan dalam meningkatkan peluang promosi dibandingkan Junior.
Variabel work_hours belum signifikan, menunjukkan jumlah jam kerja tidak cukup kuat sebagai prediktor.

Kesimpulan Praktis

Promosi sangat dipengaruhi oleh level pekerjaan dan departemen kerja.
Individu di level Senior dan di departemen Sales memiliki peluang promosi yang jauh lebih besar.
Model ini cukup baik untuk memprediksi peluang promosi, namun bisa ditingkatkan dengan menambahkan variabel lain seperti masa kerja, performa, atau pelatihan.
Treat sebagai Rasio (Numeric Rank)

sim_data_numeric <- sim_data %>%
  mutate(
    job_level_numeric = case_when(
      job_level == "Junior" ~ 1,
      job_level == "Mid" ~ 2,
      job_level == "Senior" ~ 3,
    )
)
model_numeric <- glm(promoted ~ department + job_level_numeric + work_hours, data = sim_data_numeric, family = binomial)
summary(model_numeric)

## 
## Call:
## glm(formula = promoted ~ department + job_level_numeric + work_hours, 
##     family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -2.16038    1.00967  -2.140  0.03238 *  
## departmentIT       0.76091    0.27820   2.735  0.00624 ** 
## departmentSales    1.56754    0.34885   4.493 7.01e-06 ***
## job_level_numeric  1.29532    0.23205   5.582 2.38e-08 ***
## work_hours         0.02745    0.02363   1.162  0.24528    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 477.40  on 499  degrees of freedom
## Residual deviance: 407.81  on 495  degrees of freedom
## AIC: 417.81
## 
## Number of Fisher Scoring iterations: 5

Model menggunakan job level sebagai angka bertingkat 1-3. Setiap kenaikan satu tingkat jabatan (job_level_numeric) meningkatkan peluang promosi sekitar 3,65 (exp(1.29532)) kali lipat untuk setiap level jabatan yang lebih tinggi.

Perbandingan Model

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

## $AIC_Nominal
## [1] 419.5276
## 
## $AIC_Numeric
## [1] 417.8054

Berdasarkan hasil di atas, terlihat bahwa nilai AIC untuk model dengan job level sebagai numeric (417.81) lebih rendah dari yang nominal (419.53), maka perlakuan job level sebagai numeric lebih cocok untuk data ini.

Goodness-of-Fit Model

nullmod <- glm(promoted ~ 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.1463692 (df=6)
## 
## $McFadden_R2_Numeric
## 'log Lik.' 0.1457872 (df=5)

McFadden’s \(R^2\) pada kedua model hanya berbeda sedikit sehingga tidak dapat langsung disimpulkan model yang lebih baik dalam menjelaskan data.

Visualisasi Prediksi

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 = job_level)) +
  geom_point(alpha = 0.6) +
  labs(title = "Prediksi Probabilitas (Ordinal sebagai Nominal)", x = "Work Hours", y = "Prediksi Probabilitas")

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

Visualisasi hubungan antara jam kerja (work hours) dan probabilitas mendapatkan promosi berdasarkan tingkat pekerjaan (job level) dengan dua pendekatan perlakuan ordinal cukup mirip satu sama lain. Artinya, dalam kasus ini, perlakuan job level sebagai nominal maupun numerik menghasilkan pola prediksi yang serupa, sehingga tidak terlalu memengaruhi hasil model secara drastis.

Ringkasan Koefisien Model

# Ringkasan model nominal
library(broom)
library(tidyverse)

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

library(kableExtra)
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
library(broom)
tidy(model_numeric) %>%
  kable(format = "latex", booktabs = TRUE, caption = "Ringkasan Koefisien Model dengan Ordinal sebagai Numeric") %>%
  kable_styling(latex_options = c("hold_position", "striped"))

Model Nominal: variabel department dan job level secara signifikan memengaruhi kemungkinan karyawan dipromosikan, sementara jam kerja (work_hours) tidak signifikan secara statistik (p > 0.05).
Model Numerik: variabel job level sangat signifikan berpengaruh terhadap peluang untuk dipromosikan, department juga berpengaruh, sementara jam kerja (work_hours) tetap tidak signifikan.

VIII. Pemilihan Model Regresi Logistik dan Evaluasi

Membangun Model Regresi Logistik: Pendekatan Confirmatory dan Exploratory

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

Confirmatory (Pendekatan Konfirmatori)

Pendekatan ini digunakan ketika peneliti telah memiliki teori atau hipotesis yang jelas mengenai efek atau hubungan antara variabel prediktor dan respon.

Ciri-ciri:

Model dibangun berdasarkan teori atau hasil penelitian sebelumnya.
Tujuan utamanya adalah menguji apakah efek tersebut benar-benar signifkan, bukan sekadar mencari model terbaik.
Peneliti biasanya menyusun model penuh terlebih dahulu, lalu menguji apakah penambahan atau pengurangan suatu efek (misalnya, interaksi) memberikan peningkatan model secara signifkan.
Uji signifkansi dilakukan dengan membandingkan model dengan efek tertentu dan model tanpa efek tersebut, misalnya dengan Likelihood Ratio Test.

Contoh penggunaan:

Misalnya, teori menyatakan bahwa faktor x1 dan x2 mempengaruhi probabilitas seseorang membeli produk. Maka model logistik dibangun langsung dengan x1 dan x2, lalu diuji apakah kontribusi x2 benar-benar signifikan.

Exploratory (Pendekatan Eksploratori)

Pendekatan ini digunakan ketika peneliti belum memiliki teori yang pasti atau ingin mengeksplorasi hubungan potensial antar variabel.

Ciri-ciri:

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

Proses seleksi dilakukan melalui:

Forward Selection: Mulai dari model kosong, satu per satu variabel dimasukkan jika signifkan.
Backward Elimination: Mulai dari model penuh, variabel yang tidak signifkan dikeluarkan.
Stepwise Selection: Gabungan dari keduanya, variabel dapat masuk dan keluar secara dinamis.

Tujuan:

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

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

Simulasi Data

Misalkan terdapat dataset simulasi yang merepresentasikan kemungkinan seseorang terkena risiko diabetes (y, dengan nilai 1 = berisiko, 0 = tidak), berdasarkan:

x1: umur seseorang

x2: status kebiasaan merokok (0 = tidak merokok, 1 = merokok)

x3: indeks massa tubuh (BMI)

# Load libraries
library(knitr)
library(dplyr)
library(ggplot2)
library(MASS)

## 
## Attaching package: 'MASS'

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

library(caret)

## Loading required package: lattice

## 
## Attaching package: 'caret'

## The following object is masked from 'package:purrr':
## 
##     lift

library(pROC)

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

## 
## Attaching package: 'pROC'

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

library(DescTools)

## 
## Attaching package: 'DescTools'

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

set.seed(123)

# Simulasi data
n <- 200
x1 <- rnorm(n, mean = 45, sd = 12)        # Umur
x2 <- rbinom(n, 1, 0.4)                   # Merokok
x3 <- rnorm(n, mean = 25, sd = 4)         # BMI

# Linear predictor
lin_pred <- -7 + 0.07 * x1 + 1.3 * x2 + 0.1 * x3

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

# Variabel respons biner
y <- rbinom(n, 1, p)

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

Pemilihan Model

Model Full

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

## 
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -10.87042    1.81869  -5.977 2.27e-09 ***
## x1            0.10576    0.01884   5.614 1.98e-08 ***
## x2            1.72369    0.38951   4.425 9.63e-06 ***
## x3            0.17095    0.04906   3.485 0.000492 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 247.64  on 199  degrees of freedom
## Residual deviance: 182.40  on 196  degrees of freedom
## AIC: 190.4
## 
## Number of Fisher Scoring iterations: 5

Metode Stepwise: Forward, Backward, dan Kedua Arah

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

Evaluasi Model: ROC dan AUC

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

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

## Setting direction: controls < cases

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

auc(roc_obj)

## Area under the curve: 0.8358

Pseudo R-Squared

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

##   CoxSnell Nagelkerke   McFadden 
##  0.2783570  0.3919997  0.2634666

Tabel Klasifkasi 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 120  31
##          1  18  31
##                                           
##                Accuracy : 0.755           
##                  95% CI : (0.6894, 0.8129)
##     No Information Rate : 0.69            
##     P-Value [Acc > NIR] : 0.02614         
##                                           
##                   Kappa : 0.3922          
##                                           
##  Mcnemar's Test P-Value : 0.08648         
##                                           
##             Sensitivity : 0.5000          
##             Specificity : 0.8696          
##          Pos Pred Value : 0.6327          
##          Neg Pred Value : 0.7947          
##              Prevalence : 0.3100          
##          Detection Rate : 0.1550          
##    Detection Prevalence : 0.2450          
##       Balanced Accuracy : 0.6848          
##                                           
##        'Positive' Class : 1               
##

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

## Sensitivity Specificity 
##   0.5000000   0.8695652

Model terbaik berdasarkan AIC adalah step_both karena semua metode seleksi menghasilkan AIC yang sama (190.3953). Nilai AUC dari kurva ROC sebesar 0.84, menunjukkan kemampuan klasifikasi model yang baik. Nilai pseudo R² (CoxSnell = 0.278, Nagelkerke = 0.392, McFadden = 0.263) menunjukkan model memiliki kecocokan sedang. Akurasi model mencapai 75.5%, dengan specificity tinggi (87%) namun sensitivity rendah (50%), artinya model lebih baik dalam mengidentifikasi seseorang yang tidak berisiko terkena diabetes. Secara keseluruhan, model cukup baik, namun perlu perhatian terhadap kemampuan mendeteksi seseorang yang berisiko terkena diabetes.

Metode Perbandingan Model dalam Regresi Logistik

Tujuan

Menyajikan cara membandingkan model regresi logistik menggunakan ukuran: - Deviance - AIC (Akaike Information Criterion) - Likelihood-Ratio serta menjelaskan prinsip Parsimony dalam pemilihan model.

Simulasi Data

Misalkan terdapat dataset simulasi yang merepresentasikan kemungkinan seseorang mengalami tekanan darah tinggi (y, dengan nilai 1 = hipertensi, 0 = tidak), berdasarkan:

x1: umur (dalam tahun)

x2: status merokok (0 = tidak merokok, 1 = merokok)

x3: indeks massa tubuh (BMI)

Tujuan analisis adalah membandingkan model regresi logistik yang memprediksi risiko hipertensi berdasarkan kombinasi variabel prediktor berbeda.

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

# Simulasi data
set.seed(123)
n <- 300
x1 <- rnorm(n, mean = 45, sd = 10)
x2 <- rbinom(n, 1, 0.4)
x3 <- rnorm(n, mean = 25, sd = 4)

# Membentuk probabilitas sesuai fungsi logit
lin_pred <- -2 + 0.05 * x1 + 0.8 * x2 + 0.1 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)

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

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

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

Likelihood-Ratio Test

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

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

Prinsip Parsimony

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

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

Rumus dan Penjelasan

Rumus AIC

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

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

Rumus Deviance

\[ D = -2 \, [log \, L(model) - log \, L(model \, saturasi)] \]

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

Rumus Likelihood-Ratio

\[ G^2 = -2 \, (log \, L_0 - log \, L_1) \]

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

Evaluasi Tabel Klasifkasi dan Akurasi Model

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

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0   0   0
##          1  22 278
##                                           
##                Accuracy : 0.9267          
##                  95% CI : (0.8911, 0.9535)
##     No Information Rate : 0.9267          
##     P-Value [Acc > NIR] : 0.5564          
##                                           
##                   Kappa : 0               
##                                           
##  Mcnemar's Test P-Value : 7.562e-06       
##                                           
##             Sensitivity : 1.0000          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.9267          
##          Neg Pred Value :    NaN          
##              Prevalence : 0.9267          
##          Detection Rate : 0.9267          
##    Detection Prevalence : 1.0000          
##       Balanced Accuracy : 0.5000          
##                                           
##        'Positive' Class : 1               
##

Sensitivitas dan Spesifisitas

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

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

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

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

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

## Sensitivity Specificity 
##           1           0

Kesimpulan

Berdasarkan deviance dan AIC, model 3 merupakan model yang paling baik dan cocok karena memiliki deviance dan AIC yang paling kecil.
Berdasarkan uji Likelihood Ratio Test yang membandingkan Model 2 dan Model 3, menunjukkan bahwa penambahan x3 (BMI) tidak signifikan (p = 0.127) sedangkan penambahan x2 (merokok) dari Model 1 ke Model 2 signifikan (p = 0.004). Artinya, walau Model 3 sedikit lebih baik, peningkatan fit-nya tidak cukup signifikan secara statistik dibanding Model 2.
Confusion matrix menunjukkan akurasi tinggi (92.7%), namun model ini memprediksi semua kasus sebagai positif hipertensi (sensitivitas 100%, spesifisitas 0%).
Model kurang mampu membedakan kasus negatif hipertensi dengan baik dan perlu evaluasi atau penyempurnaan agar prediksi menjadi lebih seimbang.

Penjelasan Kurva ROC (Receiver Operating Characteristic)

Kurva ROC adalah alat visual yang digunakan untuk mengevaluasi performa model klasifkasi biner. Kurva ini menunjukkan trade-of antara True Positive Rate (Sensitivity) dan False Positive Rate (1 - Specifcity) pada berbagai threshold klasifkasi.

Simulasi Data

Misalkan terdapat dataset simulasi yang menggambarkan kemungkinan seseorang mengalami hipertensi (dengan nilai 1 = hipertensi dan 0 = tidak) berdasarkan beberapa variabel prediktor:

x1 = umur, x2 = status merokok, x3 = BMI

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

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

## Setting direction: controls < cases

plot(roc_obj, main = "Kurva ROC Model Regresi Logistik")

auc(roc_obj)

## Area under the curve: 0.8686

Simulasi Pemilihan Threshold Optimal

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

thresholds <- seq(0.1, 0.9, by = 0.05)
results <- data.frame(Threshold = thresholds)

results$Sensitivity <- sapply(thresholds, function(t) {
  pred_class <- ifelse(pred >= t, 1, 0)
  cm <- table(Pred = pred_class, Obs = data$y)
  TP <- cm["1", "1"]
  FN <- cm["0", "1"]
  TP / (TP + FN)
})

results$Specificity <- sapply(thresholds, function(t) {
  pred_class <- ifelse(pred >= t, 1, 0)
  cm <- table(Pred = pred_class, Obs = data$y)
  TN <- cm["0", "0"]
  FP <- cm["1", "0"]
  TN / (TN + FP)
})

print(results)

##    Threshold Sensitivity Specificity
## 1       0.10  0.91489362   0.5947712
## 2       0.15  0.85106383   0.6862745
## 3       0.20  0.80851064   0.7320261
## 4       0.25  0.76595745   0.7712418
## 5       0.30  0.72340426   0.8104575
## 6       0.35  0.68085106   0.8366013
## 7       0.40  0.61702128   0.8954248
## 8       0.45  0.59574468   0.9150327
## 9       0.50  0.51063830   0.9281046
## 10      0.55  0.51063830   0.9477124
## 11      0.60  0.42553191   0.9607843
## 12      0.65  0.36170213   0.9738562
## 13      0.70  0.29787234   0.9803922
## 14      0.75  0.19148936   0.9869281
## 15      0.80  0.12765957   0.9869281
## 16      0.85  0.06382979   1.0000000
## 17      0.90  0.02127660   1.0000000

Precision-Recall Curve (PR Curve)

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

Visualisasi PR Curve

library(PRROC)

## Loading required package: rlang

## 
## Attaching package: 'rlang'

## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice

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

Pseudo R-squared pada Regresi Logistik

Tujuan

Menjelaskan dan menghitung pseudo R-squared dalam regresi logistik:

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

Simulasi Data

Misalkan terdapat dataset simulasi yang merepresentasikan kemungkinan seseorang membeli suatu produk (y = 1) atau tidak (y = 0), berdasarkan tiga variabel prediktor:

x1: skor minat pelanggan terhadap produk (variabel numerik)
x2: apakah pelanggan pernah membeli produk serupa sebelumnya (biner: 0 = belum, 1 = pernah)
x3: jumlah rata-rata pengeluaran pelanggan dalam sebulan (variabel numerik)

set.seed(123)
n <- 300
x1 <- rnorm(n, mean = 0, sd = 1)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n, mean = 0, sd = 1)

lin_pred <- -2 + 1.5 * x1 + 0.9 * x2 + 0.4 * x3
p <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, p)
data <- data.frame(y = as.factor(y), x1, x2, x3)

Model Logistik dan Null Model

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

Likelihood dan Rumus

\[ R^2_{Cox\;and\;Snell} = 1-(\frac{L_0}{L_M})^{2/n} \\ R^2_{McFadden} = 1-\frac{log\,L_M}{log\,L_0} \]

Dengan:

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

Perhitungan Manual R-squared

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

Perhitungan Otomatis dengan Package Tambahan

Menggunakan pscl

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

## Loading required package: pscl

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

pR2(model)

## fitting null model for pseudo-r2

##          llh      llhNull           G2     McFadden         r2ML         r2CU 
## -135.6160076 -175.9638404   80.6956658    0.2292962    0.2358457    0.3415127

Menggunakan rcompanion

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

## Loading required package: rcompanion

nagelkerke(model)

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

Menggunakan DescTools

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

##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.2292962       0.2065642       0.2358457       0.3415127       0.2119690 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.3926615       0.2549929       0.3785301       0.2593478     279.2320151 
##             BIC          logLik         logLik0              G2 
##     294.0471450    -135.6160076    -175.9638404      80.6956658

Interpretasi

McFadden R² = 0.229 → Lebih dari 0.2, artinya model memiliki kecocokan yang baik.
Cox & Snell R² = 0.236 → Nilai ini konservatif dan tidak bisa mencapai 1, namun tetap menunjukkan kemampuan prediktif yang wajar.
Nagelkerke R² = 0.342 → Sudah cukup tinggi, menandakan model cukup baik dalam menjelaskan variabilitas data.

Ketiganya menunjukkan bahwa model memiliki kekuatan prediktif yang layak dan cocok digunakan.

IX. Distribusi Multinomial

Distribusi multinomial adalah perluasan dari distribusi binomial untuk lebih dari dua kategori. Jika \(X_1, X_2, ..., X_k\) menyatakan banyaknya kejadian dalam masing-masing dari \(k\) kategori, maka:

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

dengan \(\sum_{i=1}^kx_i=n\) dan \(\sum_{i=1}^kp_i=1\).

Studi Kasus

Sebuah perusahaan melakukan survei terhadap 12 pelanggan untuk mengetahui preferensi minuman favorit dari tiga pilihan: Kopi (K), Teh (T), dan Jus (J).

Hasil survei:

Kopi: 5 orang
Teh: 4 orang
Jus: 3 orang

Probabilitas teoretik preferensi:

\(p_K\) = 0.4,
\(p_T\) = 0.3,
\(p_J\) = 0.3

Pertanyaan: Berapa peluang bahwa dari 12 orang akan ada 5 yang memilih Kopi, 4 memilih Teh, dan 3 memilih Jus?

Rumus Distribusi Multinomial

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

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

Perhitungan Manual di R

n <- 12
x <- c(5, 4, 3)
p <- c(0.4, 0.3, 0.3)

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

# Hitung peluang
peluang <- koefisien * prod(p^x)
peluang

## [1] 0.06207861

Probabilitas bahwa 5 orang memilih Kopi, 4 orang memilih Teh, dan 3 orang memilih Jus (dengan proporsi 0.4, 0.3, dan 0.3) adalah 0.06207861. Distribusi multinomial digunakan untuk menghitung peluang dari percobaan dengan hasil beberapa kategori, dan merupakan perluasan dari distribusi binomial.

Multinomial Logistic Regression

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

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

\[ log\;(\frac{P(Y=j)}{P(Y=K)})=\beta_{j_0}\,+\,\beta_{j_1x_1}\,+\,...\,+\,\beta_{j_px_p} \]

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

Contoh Kasus

Sebuah perusahaan ingin memahami faktor-faktor yang memengaruhi preferensi jenis metode kerja karyawan: WFO (Work From Office), Hybrid, atau WFH (Work From Home).

Perusahaan melakukan survei terhadap 150 karyawan dan mengumpulkan data sebagai berikut:

WorkMode: Metode kerja yang dipilih (WFO, Hybrid, WFH)
Age: Usia karyawan
Department: Divisi tempat bekerja (IT, Marketing, HR)
JobSatisfaction: Tingkat kepuasan kerja (skala 1–10)

Tujuannya adalah:

Mengetahui bagaimana usia, divisi, dan kepuasan kerja memengaruhi preferensi metode kerja karyawan.

Simulasi Data

set.seed(123)
n <- 150

Department <- sample(c("IT", "Marketing", "HR"), n, replace = TRUE)
Age <- round(rnorm(n, mean = 35, sd = 7))
JobSatisfaction <- round(pmax(pmin(rnorm(n, mean = 6, sd = 2), 10), 1))  # dibatasi antara 1 dan 10

# Simulasikan WorkMode berdasarkan probabilitas berbeda per Department
WorkMode <- sapply(Department, function(dep) {
  if (dep == "IT") {
    sample(c("WFO", "Hybrid", "WFH"), size = 1, prob = c(0.3, 0.5, 0.2))
  } else if (dep == "Marketing") {
    sample(c("WFO", "Hybrid", "WFH"), size = 1, prob = c(0.2, 0.4, 0.4))
  } else {
    sample(c("WFO", "Hybrid", "WFH"), size = 1, prob = c(0.5, 0.3, 0.2))
  }
})

df <- data.frame(WorkMode = factor(WorkMode), Age, Department = factor(Department), JobSatisfaction)
df$WorkMode <- relevel(df$WorkMode, ref = "WFO")  # baseline
head(df)

Estimasi Model

library(nnet)
model_mnlogit <- multinom(WorkMode ~ Age + Department + JobSatisfaction, data = df)

## # weights:  18 (10 variable)
## initial  value 164.791843 
## iter  10 value 153.011885
## final  value 152.930674 
## converged

summary(model_mnlogit)

## Call:
## multinom(formula = WorkMode ~ Age + Department + JobSatisfaction, 
##     data = df)
## 
## Coefficients:
##        (Intercept)         Age DepartmentIT DepartmentMarketing JobSatisfaction
## Hybrid   0.7681111 -0.02576204   0.52207900           0.8888519      0.02760348
## WFH      0.4613555 -0.02192056   0.07981353           1.6585998     -0.07627654
## 
## Std. Errors:
##        (Intercept)        Age DepartmentIT DepartmentMarketing JobSatisfaction
## Hybrid    1.175720 0.02723358    0.4734180           0.4990938      0.09921931
## WFH       1.356672 0.03209973    0.6147293           0.5549002      0.11676194
## 
## Residual Deviance: 305.8613 
## AIC: 325.8613

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 DepartmentIT DepartmentMarketing JobSatisfaction
## Hybrid      0.5136 0.3442       0.2701              0.0749          0.7809
## WFH         0.7338 0.4947       0.8967              0.0028          0.5136

Interpretasi:

Koefisien untuk kategori “Hybrid” dan “WFH” dibandingkan dengan baseline “WFO”.
Variabel DepartmentMarketing berpengaruh signifikan terhadap preferensi WFH (p-value = 0.0028 < 0.05).
Variabel lain tidak signifikan karena p-value > 0.05.

Prediksi dan Validasi

df$Predicted <- predict(model_mnlogit)
table(Predicted = df$Predicted, Actual = df$WorkMode)

##          Actual
## Predicted WFO Hybrid WFH
##    WFO     16      9   4
##    Hybrid  23     48  25
##    WFH      5     10  10

Kesimpulan

Model regresi logistik multinomial berhasil digunakan untuk:

Menganalisis hubungan antara atribut karyawan dan preferensi sistem kerja (WFO, Hybrid, WFH)
Mengetahui faktor signifikan yang memengaruhi pilihan sistem kerja
Memungkinkan prediksi sistem kerja yang dipilih oleh karyawan baru berdasarkan karakteristiknya

X. Regresi Logistik Ordinal

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

Model ini berbeda dengan:

Regresi logistik biner: hanya 2 kategori
Regresi logistik multinomial: kategori > 2 tetapi tidak berurutan

Konsep Cumulative Logit Model

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

\[ log\;(\frac{P(Y\leq{j})}{P(Y>j)})\;=\;\alpha_j\;+\;\beta_x \]

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

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

Interpretasi Koefisien

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

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

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

Odds ratio:

\[ OR=e^\beta \]

Contoh Kasus

Misalkan kita memiliki data fiktif tingkat kepedulian seseorang terhadap lingkungan (1: Rendah, 2: Sedang, 3: Tinggi).

set.seed(123)
n <- 200
activity_freq <- round(runif(n, 0, 10))  # frekuensi kegiatan 0–10 kali/bulan
score <- 4 + 0.6 * activity_freq + rnorm(n)  # skor kepedulian (nilai acak)

care_level <- cut(score,
                  breaks = c(-Inf, 5.5, 7.5, Inf),
                  labels = c("Rendah", "Sedang", "Tinggi"),
                  ordered_result = TRUE)

df <- data.frame(care_level, activity_freq)
head(df)

Estimasi Model Ordinal

model_ord <- polr(care_level ~ activity_freq, data = df, Hess = TRUE)
summary(model_ord)

## Call:
## polr(formula = care_level ~ activity_freq, data = df, Hess = TRUE)
## 
## Coefficients:
##               Value Std. Error t value
## activity_freq 1.039     0.1086   9.561
## 
## Intercepts:
##               Value  Std. Error t value
## Rendah|Sedang 2.4519 0.3938     6.2254 
## Sedang|Tinggi 6.1368 0.6522     9.4097 
## 
## Residual Deviance: 235.4224 
## AIC: 241.4224

Nilai P-Value

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

##                  Value Std. Error  t value
## activity_freq 1.038721  0.1086469 9.560515
## Rendah|Sedang 2.451874  0.3938479 6.225435
## Sedang|Tinggi 6.136831  0.6521807 9.409709

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
## activity_freq 1.038721  0.1086469 9.560515       0
## Rendah|Sedang 2.451874  0.3938479 6.225435       0
## Sedang|Tinggi 6.136831  0.6521807 9.409709       0

Prediksi Probabilitas

newdata <- data.frame(activity_freq = 3:7)
predict(model_ord, newdata = newdata, type = "probs")

##        Rendah    Sedang    Tinggi
## 1 0.339777082 0.6137221 0.0465008
## 2 0.154072747 0.7248160 0.1211113
## 3 0.060555456 0.6591957 0.2802488
## 4 0.022303625 0.4538412 0.5238552
## 5 0.008008821 0.2353755 0.7566157

Interpretasi

Koefisien activity_freq = 1.039: Artinya, semakin sering seseorang ikut kegiatan lingkungan, peluangnya untuk berada di kategori kepedulian lebih tinggi (Sedang atau Tinggi) semakin besar.
Semua p-value = 0 → Sangat signifikan. Variabel activity_freq memang berpengaruh terhadap tingkat kepedulian.
Saat activity_freq = 3: probabilitas “Rendah” = 34%, “Tinggi” = 4%
Saat activity_freq = 7: probabilitas “Tinggi” naik jadi 76%
Semakin tinggi aktivitas, semakin besar peluang tingkat kepedulian dikategorikan sebagai “Tinggi”.

XI. Log Linear Model

Tabel Kontingensi

Tabel kontingensi menyajikan jumlah frekuensi dari kombinasi kategori antar variabel.

Contoh tabel 2x2:

table_data <- matrix(c(45, 25, 15, 65), nrow = 2,
                     dimnames = list(Merokok = c("Ya", "Tidak"),
                                     Penyakit = c("Ya", "Tidak")))

table_data

##        Penyakit
## Merokok Ya Tidak
##   Ya    45    15
##   Tidak 25    65

Model Loglinear

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

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

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

library(MASS)
loglm(~ Merokok * Penyakit, data = table_data)

## Call:
## loglm(formula = ~Merokok * Penyakit, 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\,(\frac{p}{1-p})=\beta_0\,+\,\beta_1x \]

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

Contoh:

data_glm <- data.frame(
Penyakit = c(1, 0, 1, 0),
Merokok = factor(c("Ya", "Ya", "Tidak", "Tidak")),
Frek = c(45, 25, 15, 65)
)
model_logit <- glm(Penyakit ~ Merokok, weights = Frek, family = binomial, data = data_glm)
summary(model_logit)

## 
## Call:
## glm(formula = Penyakit ~ Merokok, family = binomial, data = data_glm, 
##     weights = Frek)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.4663     0.2864  -5.119 3.07e-07 ***
## MerokokYa     2.0541     0.3798   5.408 6.37e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 201.90  on 3  degrees of freedom
## Residual deviance: 168.46  on 2  degrees of freedom
## AIC: 172.46
## 
## Number of Fisher Scoring iterations: 5

Tabel Kontingensi dan Model Loglinier

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

data <- matrix(c(45, 25, 15, 65), nrow = 2,
               dimnames = list(Merokok = c("Ya", "Tidak"),
                               Penyakit = c("Ya", "Tidak")))
data

##        Penyakit
## Merokok Ya Tidak
##   Ya    45    15
##   Tidak 25    65

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

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

Model Saturated

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

Cocok sempurna terhadap data
Tidak mengasumsikan independensi antar variabel

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

model_saturated <- loglm(~ Merokok * Penyakit, data = data)
summary(model_saturated)

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

Model Independent

Model independen mengasumsikan bahwa tidak ada interaksi antara variabel:

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

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

model_indep <- loglm(~ Merokok + Penyakit, data = data)
summary(model_indep)

## Formula:
## ~Merokok + Penyakit
## attr(,"variables")
## list(Merokok, Penyakit)
## attr(,"factors")
##          Merokok Penyakit
## Merokok        1        0
## Penyakit       0        1
## attr(,"term.labels")
## [1] "Merokok"  "Penyakit"
## 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 33.44517  1 7.330155e-09
## Pearson          32.25446  1 1.352460e-08

Odds Ratio dan Interpretasi

Odds ratio untuk tabel 2x2:

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

Interpretasi nilai OR:

OR = 1: Tidak ada asosiasi
OR > 1: Asosiasi positif
OR < 1: Asosiasi negatif

Estimasi Parameter

Dalam model saturated:

Estimasi dilakukan dengan pembatasan seperti sum-to-zero
Estimasi parameter dilakukan dengan iterative proportional ftting (IPF)

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

## [1] 2.054124

Model Lebih Sederhana dan Perbandingan Model

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

anova(model_indep, model_saturated)

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

Studi Kasus

Sebuah survei hipotetik dilakukan untuk melihat hubungan antara frekuensi konsumsi sayur dan tingkat energi harian seseorang.

# Buat tabel kontingensi 3x2
data_survey <- matrix(c(25, 85,
                        60, 140,
                        80, 210),
                      nrow = 3, byrow = TRUE,
                      dimnames = list(Konsumsi = c("Jarang", "Cukup", "Sering"),
                                      Energi = c("Rendah", "Tinggi")))

# Tampilkan tabel
ftable(data_survey)

##          Energi Rendah Tinggi
## Konsumsi                     
## Jarang              25     85
## Cukup               60    140
## Sering              80    210

loglm(~ Konsumsi + Energi, data = data_survey)

## Call:
## loglm(formula = ~Konsumsi + Energi, data = data_survey)
## 
## Statistics:
##                       X^2 df  P(> X^2)
## Likelihood Ratio 1.924125  2 0.3821039
## Pearson          1.884809  2 0.3896898

Interpretasi:

Nilai p-value dari Likelihood Ratio dan Pearson masing-masing adalah 0.3821 dan 0.3897, yang lebih besar dari 0.05. Artinya, tidak ada bukti yang signifikan untuk menyatakan bahwa ada hubungan antara Konsumsi Sayur dan Tingkat Energi sehingga variabel konsumsi sayur dan energi harian bersifat independen dalam data ini.

XII. Referensi

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

Azis, I., Sumertajaya, I. M., Purwaningsih, S. S., & Tjahjawati, S. S. (2023). Penentuan faktor kemiskinan Indonesia menggunakan regresi logistik ordinal. Journal of Mathematics, Computations, and Statistics, 6(1), 61‑65.

Cahyadhie, A. (n.d.). Buku Analisis Regresi Logistik [PDF].

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

Harlan, J. (2018). Analisis regresi logistik. Gunadarma.

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

Maryana, M. (2013). Model log linier yang terbaik untuk analisis data kualitatif pada tabel kontingensi tiga arah. Malikussaleh Industrial Engineering Journal, 2(2), 3–10.

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

Nurrizqi, A. I., Erfiani, I., Indahwati, & Fitrianto, A. (2022). Pemodelan regresi logistik berbasis backward elimination untuk mengetahui faktor yang memengaruhi tingkat kemiskinan di Indonesia tahun 2021. Jurnal Statistika dan Aplikasinya, 6(2), 160‑170.

Rahmawati. (–). Regresi logistik untuk kasus anemia ibu hamil. Dalam buku Statistik Kesehatan.

Roflin, E., Riana, F., Munarsih, E., Pariyana, & Andriyani, I. (2023). Regresi logistik biner dan multinomial. Penerbit NEM.

Tirta, I. M. (2015). Theori model linier tergeneralisir (GLM) dengan variabel kualitatif (Dummy), natural spline dan B-spline. Laboratorium Statistika, FMIPA Universitas Jember.

E-Book Analisis Data Kategori

Gina Kustiana

2025-04-12

KATA PENGANTAR

I. Tabel Kontingensi 2 x 2

Contoh Kasus

Distribusi Peluang

- Peluang Bersama (Joint Probability)

- Peluang Marginal

- Peluang Bersyarat (Conditional Probability)

Ukuran Asosiasi

- Beda Peluang (Risk Difference)

- Risiko Relatif (Relative Risk)

- Rasio Odds (Odds Ratio)

II. Inferensi Tabel Kontingensi Dua Arah

Estimasi

- Estimasi Titik

- Estimasi Interval

Uji Hipotesis

- Uji Proporsi

- Uji Asosiasi

- Uji Independensi

Analisis Residual dalam Tabel Kontingensi

- Analisis Pearson Residual

- Deteksi Outlier Menggunakan Residual

III. Tabel Kontingensi Tiga Arah

Contoh Kasus

Tabel Parsial & Tabel Marginal

- Tabel Parsial

- Tabel Marginal

Distribusi Peluang

- Peluang Bersama (Joint Probability)

- Peluang Marginal

- Peluang Bersyarat (Conditional Probability)

Tabel Peluang Bersyarat

Ukuran Asosiasi

- Beda Peluang (Risk Difference)

- Risiko Relatif (Relative Risk)

- Rasio Odds (Odds Ratio)

Conditional Independence

IV. Inferensi Tabel Kontingensi Tiga Arah

Pengujian Independensi

Metode Cochran-Mantel-Haenszel (CMH)

V. Generalized Linear Model (GLM)

Exponential Family (Keluarga Eksponensial)

Contoh

Model Regresi Logistik

Contoh

Model Regresi Poisson

Contoh

VI. Inferensi Generalized Linear Model (GLM)

Inferensi Parameter

- Uji Wald

- Uji Likelihood Ratio (Chi-Square)

VII. Regresi Logistik dengan Prediktor Nominal, Ordinal, dan Rasio

Contoh Kasus

VIII. Pemilihan Model Regresi Logistik dan Evaluasi

Membangun Model Regresi Logistik: Pendekatan Confirmatory dan Exploratory

Metode Stepwise: Forward, Backward, dan Kedua Arah

Evaluasi Model: ROC dan AUC

Pseudo R-Squared

Tabel Klasifkasi dan Evaluasi

Metode Perbandingan Model dalam Regresi Logistik

Likelihood-Ratio Test

Prinsip Parsimony

Evaluasi Tabel Klasifkasi dan Akurasi Model

Penjelasan Kurva ROC (Receiver Operating Characteristic)

Precision-Recall Curve (PR Curve)

Pseudo R-squared pada Regresi Logistik

IX. Distribusi Multinomial

Studi Kasus

Multinomial Logistic Regression

Contoh Kasus

X. Regresi Logistik Ordinal

Konsep Cumulative Logit Model

Interpretasi Koefisien

Contoh Kasus

XI. Log Linear Model

Tabel Kontingensi dan Model Loglinier

Model Saturated