1. Analisis Data Kategori

1.1 Pendahuluan

Dalam bidang statistik, data diklasifikasikan menjadi beberapa jenis berdasarkan karakteristiknya, salah satunya adalah data kategorik. Data kategorik merupakan data yang menyatakan atribut atau kategori tertentu, bukan nilai numerik yang dapat diukur secara langsung. Contoh dari data ini meliputi jenis kelamin (laki-laki atau perempuan), status pernikahan (menikah, belum menikah, atau bercerai), serta tingkat pendidikan (SMA, Diploma, Sarjana, dan sebagainya). Karena sifatnya yang bersifat kualitatif, data kategorik memerlukan metode analisis yang berbeda dibandingkan dengan data numerik.

Analisis data kategorik bertujuan untuk menggali informasi dari data yang bersifat kualitatif, dengan cara mengeksplorasi distribusi, hubungan, serta pengaruh antar variabel. Beberapa teknik yang umum digunakan dalam analisis ini meliputi penyusunan tabel kontingensi, pengujian independensi antar variabel dengan uji chi-square, serta pemodelan dengan regresi logistik untuk variabel dependen kategorik. Pendekatan ini sangat bermanfaat dalam berbagai disiplin ilmu, seperti ilmu sosial, kesehatan, pemasaran, dan pendidikan, di mana data kategorik sering kali menjadi jenis data utama.

1.2 Tujuan

Mengetahui distribusi dan proporsi kategori
Menyajikan sebaran data dalam setiap kategori untuk memahami komposisi atau frekuensi kemunculan masing-masing kategori.

Mengidentifikasi hubungan antar variabel kategorik
Menganalisis apakah ada keterkaitan atau asosiasi antara dua atau lebih variabel kategorik.

Menguji signifikansi hubungan antar kategori
Menggunakan uji hipotesis seperti uji chi-square untuk menguji apakah hubungan antar kategori signifikan secara statistik.

Membangun model prediksi untuk variabel kategorik
Menggunakan model seperti regresi logistik untuk memprediksi kategori berdasarkan variabel lain.

Memberikan pemahaman lebih dalam terhadap fenomena kualitatif
Menganalisis data untuk memberikan wawasan yang lebih mendalam tentang fenomena berdasarkan data kategorik.

2. Metode dalam Analisis Data Kategori

2.1 Tabel Kontingensi dan Uji Chi Square

Tujuan: Untuk menggambarkan distribusi frekuensi dari dua atau lebih variabel kategorik dan mengidentifikasi pola hubungan di antara variabel tersebut. Tabel ini sering digunakan untuk analisis hubungan antara variabel, seperti jenis kelamin dan status pekerjaan.

# Membuat data tabel kontigensi
data <- data.frame(
  Gender = c("Male", "Female", "Female", "Male", "Female", "Male", "Female", "Male"),
  JobStatus = c("Employed", "Employed", "Unemployed", "Employed", "Unemployed", "Employed", "Unemployed", "Employed")
)

# Membuat tabel kontingensi
table_data <- table(data$Gender, data$JobStatus)
print(table_data)

##         
##          Employed Unemployed
##   Female        1          3
##   Male          4          0

Uji Chi-Square (Chi-Square Test)
Tujuan: Untuk menguji apakah terdapat hubungan yang signifikan secara statistik antara dua variabel kategorik. Uji chi-square menguji apakah distribusi frekuensi yang diamati berbeda dengan distribusi yang diharapkan berdasarkan hipotesis nol (misalnya, tidak ada hubungan antara dua variabel).

# Membuat data baru untuk uji chi-square
data <- data.frame(
  Gender = c("Male", "Female", "Male", "Female", "Male", "Female", 
             "Male", "Female", "Male", "Female"),
  ColorPreference = c("Blue", "Red", "Blue", "Red", "Blue", "Blue", 
                      "Green", "Red", "Green", "Green")
)

# Membuat tabel kontingensi
table_data5 <- table(data$Gender, data$ColorPreference)
print(table_data5)

##         
##          Blue Green Red
##   Female    1     1   3
##   Male      3     2   0

# Uji Chi-Square
chisq_test5 <- chisq.test(table_data5)

## Warning in chisq.test(table_data5): Chi-squared approximation may be incorrect

print(chisq_test5)

## 
##  Pearson's Chi-squared test
## 
## data:  table_data5
## X-squared = 4.3333, df = 2, p-value = 0.1146

2.2 Regresi Logistik (Logistic Regression)

Tujuan: Untuk memodelkan hubungan antara variabel dependen kategorik (misalnya, dua kategori seperti sukses/gagal) dengan satu atau lebih variabel independen. Ini digunakan untuk prediksi probabilitas suatu kejadian berdasarkan faktor-faktor lain.

# Membuat data  untuk regresi logistik
data <- data.frame(
  Outcome = c(1, 0, 1, 0, 1, 1, 0, 1, 0, 1),  # 1 = Lulus, 0 = Gagal
  StudyMethod = factor(c("Online", "Offline", "Online", "Offline", "Online", "Offline", 
                         "Online", "Offline", "Online", "Offline"))
)

# Membuat model regresi logistik
logistic_model <- glm(Outcome ~ StudyMethod, family = binomial, data = data)
summary(logistic_model)

## 
## Call:
## glm(formula = Outcome ~ StudyMethod, family = binomial, data = data)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)
## (Intercept)       4.055e-01  9.129e-01   0.444    0.657
## StudyMethodOnline 5.733e-16  1.291e+00   0.000    1.000
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 13.46  on 9  degrees of freedom
## Residual deviance: 13.46  on 8  degrees of freedom
## AIC: 17.46
## 
## Number of Fisher Scoring iterations: 4

Analisis Varians Kategorik (ANOVA untuk data kategorik)
Tujuan: Untuk membandingkan proporsi kategori di antara lebih dari dua kelompok dan menguji apakah perbedaan antar kelompok tersebut signifikan. Meskipun ANOVA lebih umum digunakan untuk data numerik, ada juga aplikasi khusus untuk data kategorik.

# Membuat data baru untuk ANOVA
data <- data.frame(
  Education = factor(c("High School", "Bachelor", "Master", "High School", 
                       "Master", "Bachelor", "High School", "Bachelor", 
                       "Master", "Master")),
  Income = c(3000, 4500, 5500, 3200, 5800, 4700, 2900, 5000, 6000, 6300)
)

# Melakukan uji ANOVA
anova_result3 <- aov(Income ~ Education, data = data)
summary(anova_result3)

##             Df   Sum Sq Mean Sq F value   Pr(>F)    
## Education    2 14095667 7047833   96.11 8.13e-06 ***
## Residuals    7   513333   73333                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Koefisien Asosiasi (Cramér’s V, Phi Coefficient)
Tujuan: Untuk mengukur kekuatan asosiasi antara dua variabel kategorik. Koefisien ini memberikan ukuran besarnya hubungan yang ada, dengan nilai mendekati 0 menunjukkan hubungan yang sangat lemah dan nilai mendekati 1 menunjukkan hubungan yang sangat kuat.

library(vcd)

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

## Loading required package: grid

# Membuat data baru untuk koefisien asosiasi
data <- data.frame(
  Gender = c("Male", "Female", "Male", "Female", "Male", "Female", "Male", 
             "Female", "Female", "Male"),
  MoviePreference = c("Action", "Romance", "Action", "Romance", "Action", 
                      "Romance", "Action", "Romance", "Action", "Romance")
)

# Membuat tabel kontingensi
table_data4 <- table(data$Gender, data$MoviePreference)

# Menghitung koefisien Cramér's V
cramer_v <- assocstats(table_data4)
cramer_v$cramer

## [1] 0.6

3. Distribusi Probabilitas dalam Data Kategori

3.1 Distribusi Bernoulli

Distribusi Bernoulli adalah distribusi probabilitas untuk sebuah variabel acak yang hanya memiliki dua hasil, yaitu sukses (1) atau gagal (0). Distribusi ini digunakan untuk model yang hanya memiliki dua kategori.

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

Di mana:

$p$ adalah probabilitas terjadinya sukses (misalnya, $p = P(X = 1)$).
$1 - p$ adalah probabilitas terjadinya kegagalan (misalnya, $1 - p = P(X = 0)$).

# Probabilitas sukses (misalnya, probabilitas lulus ujian)
p <- 0.75

# Menghitung probabilitas untuk sukses (X = 1) dan gagal (X = 0)
prob_success <- p
prob_failure <- 1 - p

# Menampilkan hasil
cat("P(X = 1) =", prob_success, "\n")

## P(X = 1) = 0.75

cat("P(X = 0) =", prob_failure, "\n")

## P(X = 0) = 0.25

Penjelasan:

Dalam contoh ini, p = 0.75 yang berarti ada 75% peluang sukses (misalnya lulus ujian). Probabilitas kegagalan (X = 0) adalah 1 - 0.75 = 0.25.

3.2 Distribusi Multinomial

Distribusi Multinomial adalah generalisasi dari distribusi Bernoulli untuk variabel acak yang dapat mengambil lebih dari dua kategori. Distribusi ini sering digunakan untuk menghitung probabilitas dalam eksperimen yang memiliki lebih dari dua hasil, seperti melempar dadu yang memiliki lebih dari dua sisi.

Rumus Distribusi Multinomial:

Untuk $n$ percobaan dan $k$ kategori, distribusi multinomial memberikan probabilitas
$P(X_1 = x_1, X_2 = x_2, ..., X_k = x_k)$ sebagai:

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

Di mana:

$n$ adalah jumlah total percobaan.
$x_i$ adalah jumlah hasil untuk kategori ke-$i$.
$p_i$ adalah probabilitas kategori ke-$i$.

# Probabilitas untuk tiap kategori
p <- c(0.2, 0.3, 0.5)

# Jumlah percobaan (misalnya, melempar dadu 10 kali)
n <- 10

# Menghitung probabilitas distribusi multinomial
multinom_prob <- dmultinom(x = c(2, 3, 5), prob = p)
cat("Probabilitas untuk distribusi multinomial:", multinom_prob, "\n")

## Probabilitas untuk distribusi multinomial: 0.08505

Penjelasan:

Di sini, kita menguji distribusi probabilitas ketika kita melempar “dadu” tiga kali, dengan probabilitas untuk masing-masing kategori (misalnya, angka 1, 2, dan 3) masing-masing 0.2,0.3, dan 0.5
dmultinom() digunakan untuk menghitung probabilitas dalam distribusi multinomial berdasarkan hasil yang diberikan (misalnya, 2 angka pertama, 3 angka kedua, dan 5 angka ketiga).

3.3 Distribusi Binomial

Distribusi Binomial adalah distribusi untuk variabel acak yang dapat menghasilkan dua hasil (sukses atau gagal) dalam sejumlah percobaan tetap. Distribusi ini sangat berguna untuk mengukur jumlah keberhasilan dalam n percobaan dengan probabilitas sukses p.

Rumus Distribusi Binomial:

Distribusi Binomial memberikan probabilitas untuk mendapatkan x sukses dalam n percobaan:

Distribusi Binomial memberikan probabilitas untuk mendapatkan $x$ sukses dalam $n$ percobaan:

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

Di mana:

$\binom{n}{x}$ adalah kombinasi dari $n$ percobaan yang diambil $x$ sukses.
$p$ adalah probabilitas sukses.

# Probabilitas sukses
p <- 0.5

# Jumlah percobaan
n <- 10

# Menghitung probabilitas mendapatkan 3 sukses dalam 10 percobaan
prob_3_success <- dbinom(3, size = n, prob = p)
cat("P(X = 3) =", prob_3_success, "\n")

## P(X = 3) = 0.1171875

Penjelasan:

Di sini, kita menguji distribusi binomial dengan probabilitas sukses 50% (p = 0.5) dalam 10 percobaan. Probabilitas mendapatkan 3 sukses dihitung menggunakan fungsi dbinom().

3.4. Distribusi Poisson

Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam interval waktu atau ruang tertentu, dengan asumsi kejadian tersebut terjadi secara acak dan independen. Distribusi ini sangat berguna ketika rata-rata kejadian per unit waktu/spatial diketahui.

Rumus Distribusi Poisson

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

Di mana:

$\lambda$ adalah rata-rata jumlah kejadian (rate).
$x$ adalah jumlah kejadian yang diinginkan.

# Rata-rata jumlah kejadian (misalnya, 3 kejadian per jam)
lambda <- 3

# Menghitung probabilitas terjadi 4 kejadian
prob_4_occurrences <- dpois(4, lambda)
cat("P(X = 4) =", prob_4_occurrences, "\n")

## P(X = 4) = 0.1680314

Kesimpulan

Distribusi Bernoulli digunakan untuk data kategorikal dengan dua kategori (misalnya sukses dan gagal).
Distribusi Multinomial digunakan untuk data kategorikal dengan lebih dari dua kategori (misalnya hasil lemparan dadu).
Distribusi Binomial digunakan untuk data yang menunjukkan jumlah sukses dalam percobaan berulang.
Distribusi Poisson digunakan untuk menghitung jumlah kejadian dalam interval waktu atau ruang tertentu

4. Desain Sampling dalam analisis data kategori

Dalam analisis data kategorik, pemilihan desain sampling sangat penting karena akan memengaruhi cara data dikumpulkan dan bagaimana analisis dilakukan. Dua pendekatan utama dalam desain sampling untuk data kategorik adalahtolo prospective sampling dan retrospective sampling. Keduanya umum digunakan dalam studi observasional, terutama dalam konteks medis, epidemiologi, dan ilmu sosial.

4.1. Prospective Sampling

Prospective sampling merupakan pendekatan dalam pengambilan sampel di mana partisipan dipilih terlebih dahulu berdasarkan eksposur atau karakteristik tertentu, kemudian diikuti dalam jangka waktu tertentu untuk melihat bagaimana hasil atau outcome berkembang. Teknik ini banyak digunakan dalam studi kausal dan longitudinal karena memungkinkan pengamatan terhadap hubungan sebab-akibat secara lebih langsung.

4.1.1. Studi Eksperimental

Dalam konteks eksperimental, peserta dikelompokkan ke dalam perlakuan dan kontrol, biasanya melalui alokasi acak. Beberapa teknik pengambilan sampel yang sering digunakan meliputi:

Simple Random Sampling (SRS): Setiap anggota populasi memiliki peluang yang sama untuk dipilih.
Stratified Random Sampling: Populasi dibagi ke dalam strata berdasarkan atribut tertentu, lalu sampel dipilih secara acak dari tiap strata.
Cluster Sampling: Populasi dikelompokkan ke dalam unit-unit (cluster), kemudian beberapa cluster dipilih secara acak dan semua anggota dalam cluster tersebut diteliti.

4.1.2. Studi Kohort

Studi kohort merupakan bentuk observasi di mana sekelompok individu dengan karakteristik yang sama (misalnya status paparan) diikuti dalam waktu tertentu untuk mencatat kejadian yang diamati. Beberapa teknik sampling yang digunakan antara lain:

Census Sampling: Semua individu dalam populasi dimasukkan dalam studi.
Systematic Sampling: Sampel dipilih menggunakan pola atau interval tertentu dari daftar populasi.
Matched Sampling: Setiap anggota kelompok dibandingkan atau dipasangkan dengan individu serupa dari kelompok lain, berdasarkan variabel tertentu (misalnya usia, jenis kelamin).

4.2. Retrospective Sampling

Retrospective sampling adalah pendekatan di mana data dikumpulkan setelah outcome atau hasil telah terjadi. Pendekatan ini umum digunakan dalam studi observasional, terutama ketika peneliti ingin menelusuri faktor-faktor risiko yang mungkin berkontribusi terhadap suatu kejadian.

4.2.1. Studi Kasus-Kontrol

Dalam studi kasus-kontrol, sampel terdiri dari individu yang telah mengalami suatu kejadian (kasus) dan individu yang belum mengalaminya (kontrol). Beberapa teknik sampling yang lazim digunakan meliputi:

Purposive Sampling: Sampel dipilih secara sengaja karena memiliki karakteristik penting yang relevan dengan tujuan studi.
Snowball Sampling: Sampel awal diminta merekomendasikan individu lain yang memiliki karakteristik serupa.
Incidence Density Sampling: Kontrol dipilih dari populasi yang sama dengan mempertimbangkan waktu kejadian kasus, untuk menghasilkan estimasi yang lebih akurat terhadap risiko relatif.

4.2.2. Studi Kohort Retrospektif

Studi kohort retrospektif adalah jenis studi observasional di mana kelompok individu dibentuk berdasarkan data historis mengenai paparan atau karakteristik tertentu, lalu peneliti menelusuri kembali catatan atau informasi yang sudah ada untuk menentukan apakah peserta mengalami outcome tertentu. Berbeda dengan studi kohort prospektif yang dilakukan ke depan, pendekatan ini bersifat menelusuri ke belakang (retrospektif), namun tetap mempertahankan struktur kohort.

Jenis sampling yang sering diterapkan dalam studi kohort retrospektif antara lain:

Record-Based Sampling: Sampel ditentukan berdasarkan data sekunder seperti rekam medis, catatan perusahaan, atau basis data administratif yang berisi informasi tentang eksposur dan outcome.
Exposure-Based Sampling: Kelompok individu dibentuk berdasarkan informasi terdahulu mengenai status paparan, lalu dilakukan penelusuran apakah mereka mengalami outcome tertentu.
Time-Defined Sampling: Sampel diambil dari individu yang tercatat dalam suatu periode waktu tertentu di masa lalu, untuk dianalisis hubungan antara paparan dan outcome berdasarkan catatan historis.

4.3 Tabel Perbandingan

5. Tabel Kontigensi 2x2

Tabel kontingensi 2x2 adalah suatu tabel yang digunakan untuk menyajikan hubungan antara dua variabel kategorik yang masing-masing memiliki dua kategori. Tabel ini membantu dalam menganalisis asosiasi atau hubungan antara dua variabel, serta menjadi dasar untuk berbagai uji statistik seperti uji chi-square atau perhitungan odds ratio. Tabel ini sering digunakan dalam studi epidemiologi, eksperimen sosial, maupun survei, terutama saat membandingkan proporsi antar dua kelompok terhadap dua kondisi atau perlakuan.

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

Di mana:

$p$ adalah probabilitas terjadinya sukses (misalnya, $p = P(X = 1)$).
$1 - p$ adalah probabilitas terjadinya kegagalan (misalnya, $1 - p = P(X = 0)$).

Tabel Kontingensi

	Tidak Ada Logo Halal	Ada Logo Halal	Total
Beli	5	0	5
Tidak Beli	44	96	140
Total	49	96	145

5.1 Distribusi Peluang dalam Tabel Kontigensi

5.1.1 Peluang Bersama (Joint Probability)

Peluang bersama merupakan peluang dua kejadian terjadi secara bersamaan.

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

Di mana:

$n_{ij}$ adalah frekuensi kejadian gabungan dari $A_i$ dan $B_j$
$n$ adalah jumlah total observasi

5.1.2 Peluang Marginal (Marginal Probability)

Peluang marginal merupakan peluang dimana terjadinya satu kejadian tanpa mempertimbangkan kejadian lainnya.

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

Di mana:

$n_{i\cdot}$ adalah jumlah total observasi pada baris ke-$i$ (untuk $A_i$)
$n_{\cdot j}$ adalah jumlah total observasi pada kolom ke-$j$ (untuk $B_j$)
$n$ adalah jumlah total observasi

5.1.3 Peluang Bersyarat (Conditional Probability)

Peluang bersyarat adalah peluang dimana suatu kejadian terjadi dengan syarat kejadian lain sudah terjadi.

Rumus Probabilitas Bersyarat

\[ P(A_i \mid B_j) = \frac{P(A_i, B_j)}{P(B_j)} = \frac{n_{ij}}{n_{\cdot j}} \]

Di mana:

$P(A_i \mid B_j)$ adalah probabilitas terjadinya $A_i$ dengan syarat $B_j$ terjadi
$n_{ij}$ adalah jumlah observasi pada kategori $A_i$ dan $B_j$
$n_{\cdot j}$ adalah jumlah total observasi untuk kategori $B_j$

# Membuat tabel kontingensi

tabel <- matrix(c(5, 44, 0, 96),
                nrow = 2, byrow = TRUE,
                dimnames = list(
                  "Keputusan Pembelian" = c("Beli", "Tidak Beli"),
                  "Logo Halal" = c("Tidak ada", "Ada")
                ))

# Total keseluruhan
total <- sum(tabel)

# Peluang bersama = frekuensi / total
peluang_bersama <- tabel / total
print("Peluang Bersama:")

## [1] "Peluang Bersama:"

print(round(peluang_bersama, 3))

##                    Logo Halal
## Keputusan Pembelian Tidak ada   Ada
##          Beli           0.034 0.303
##          Tidak Beli     0.000 0.662

# Peluang marginal baris (Keputusan Pembelian)
peluang_marginal_baris <- rowSums(tabel) / total
print("Peluang Marginal (Baris):")

## [1] "Peluang Marginal (Baris):"

print(round(peluang_marginal_baris, 3))

##       Beli Tidak Beli 
##      0.338      0.662

# Peluang marginal kolom (Logo Halal)
peluang_marginal_kolom <- colSums(tabel) / total
print("Peluang Marginal (Kolom):")

## [1] "Peluang Marginal (Kolom):"

print(round(peluang_marginal_kolom, 3))

## Tidak ada       Ada 
##     0.034     0.966

# Peluang bersyarat: P(Beli | Tidak ada logo halal)
p_beli_diberi_tidak_ada <- tabel["Beli", "Tidak ada"] / sum(tabel[, "Tidak ada"])
# Peluang bersyarat: P(Beli | Ada logo halal)
p_beli_diberi_ada <- tabel["Beli", "Ada"] / sum(tabel[, "Ada"])

print("Peluang Bersyarat:")

## [1] "Peluang Bersyarat:"

print(paste("P(Beli | Tidak ada logo halal) =", round(p_beli_diberi_tidak_ada, 3)))

## [1] "P(Beli | Tidak ada logo halal) = 1"

print(paste("P(Beli | Ada logo halal) =", round(p_beli_diberi_ada, 3)))

## [1] "P(Beli | Ada logo halal) = 0.314"

Langkah 1: Hitung Peluang Bersama

\[ P(\text{Beli, Tidak ada}) = \frac{5}{145} = 0.0345 \] \[ P(\text{Beli,Ada}) = \frac{0}{145} = 0 \] \[ P(\text{Tidak Beli, Tidak ada}) = \frac{44}{145} = 0.3034 \] \[ P(\text{Tidak Beli,Ada}) = \frac{96}{145} = 0.0345 \]

Langkah 2: Hitung Peluang Marginal

\[ P(\text{Beli}) = \frac{5}{145} = 0.0345 \]

\[ P(\text{Tidak Beli}) = \frac{140}{145} = 0.9655 \]

\[ P(\text{Tidak ada logo halal}) = \frac{49}{145} = 0.3379 \] \[ P(\text{Ada logo halal}) = \frac{96}{145} = 0.6621 \]

Langkah 3 : Hitung Peluang Bersyarat

\[ P(\text{Beli | Tidak ada}) = \frac{5}{5} = 1 \] \[ P(\text{Beli | Ada}) = \frac{44}{140} = 0.3143 \]

Interpretasi

Peluang gabungan memberikan informasi tentang kemungkinan dua kejadian terjadi bersamaan, misalnya seseorang membeli produk dan produk itu tidak memiliki logo halal.
Peluang marginal membantu kita mengetahui seberapa besar kemungkinan dari satu kejadian tanpa melihat variabel lain.
Peluang bersyarat menjelaskan bagaimana peluang satu kejadian berubah setelah kita tahu informasi tentang kejadian lainnya.

5.2 Ukuran Asosiasi dalam Data Kategori 2 x 2

Ukuran asosiasi dalam tabel kontingensi 2x2 adalah alat statistik yang digunakan untuk mengukur kekuatan dan arah hubungan antara dua variabel kategori. Tabel ini biasanya merepresentasikan dua kelompok (misalnya “terpapar” dan “tidak terpapar”) dan dua hasil (misalnya “positif” dan “negatif”). Dengan menggunakan ukuran asosiasi, kita dapat memahami apakah keberadaan suatu faktor berhubungan dengan peningkatan atau penurunan kemungkinan suatu kejadian terjadi.

Analisis asosiasi banyak diterapkan di berbagai bidang seperti epidemiologi, pemasaran, ilmu sosial, dan kedokteran. Contohnya dalam dunia medis, asosiasi digunakan untuk mengetahui apakah ada hubungan antara konsumsi makanan tertentu dengan risiko terkena penyakit jantung. Dalam pemasaran, analisis asosiasi dapat membantu menentukan apakah kehadiran label “halal” berpengaruh terhadap keputusan pembelian konsumen terhadap suatu produk.

Risk Difference (RD) adalah selisih antara probabilitas kejadian pada kelompok terpapar dan tidak terpapar, menggambarkan seberapa besar risiko bertambah/berkurang akibat paparan.

Relative Risk (RR) adalah rasio antara probabilitas kejadian pada kelompok terpapar dibandingkan dengan kelompok tidak terpapar, menggambarkan seberapa kali lebih besar risiko pada kelompok terpapar.

Odds Ratio (OR) adalah perbandingan odds (peluang kejadian dibanding tidak kejadian) antara dua kelompok, dan sering digunakan pada studi kasus-kontrol.

Uji Chi-Square digunakan untuk menguji apakah terdapat hubungan yang signifikan secara statistik antara dua variabel kategori.

Fisher’s Exact Test digunakan sebagai alternatif uji chi-square ketika ukuran sampel kecil atau ada nilai frekuensi harapan yang sangat kecil, karena lebih akurat dalam kondisi tersebut.

5.2.1 Risk Difference (RD)

Risk Difference mengukur selisih probabilitas kejadian antara dua kelompok. Nilai RD bisa positif (risiko lebih besar pada kelompok terpapar), nol (tidak ada perbedaan), atau negatif (risiko lebih besar pada kelompok tidak terpapar).

Rumus : \[ RD = P(kejadian | terpapar) - P(kejadian | tidak terpapar) \]

# Nilai dari tabel
beli_tanpa_logo <- 5
total_tanpa_logo <- 49

beli_dengan_logo <- 0
total_dengan_logo <- 96

# Hitung risiko masing-masing kelompok
risk_terpapar <- beli_tanpa_logo / total_tanpa_logo
risk_tidak_terpapar <- beli_dengan_logo / total_dengan_logo

# Hitung Risk Difference
RD <- risk_terpapar - risk_tidak_terpapar
RD

## [1] 0.1020408

Hasil dan Interpretasi

Peluang beli pada produk tanpa logo halal:

\[ RD = \frac{5}{49} - \frac{0}{96} = 0.1020 \]

Kesimpulan:
Produk tanpa logo halal memiliki risiko dibeli sebesar 10.2% lebih tinggi dibandingkan produk yang berlogo halal.

5.2.2 Relative Risk (RR)

Relative Risk menunjukkan berapa kali lebih besar (atau kecil) risiko kejadian terjadi pada kelompok terpapar dibandingkan tidak terpapar.

Rumus :

\[ RR = \frac{P(kejadian | terpapar)}{P(kejadian | tidak\ terpapar)} \]

RR <- risk_terpapar / risk_tidak_terpapar
RR

## [1] Inf

Hasil dan interpretasi

\[ RR = \frac{0.1020}{0} \Rightarrow tidak terdefinisi (∞) \]

Karena tidak ada yang membeli di kelompok “ada logo halal”, maka Relative Risk tidak bisa dihitung dan dianggap tak terhingga. Ini menunjukkan bahwa seluruh kejadian hanya terjadi di kelompok terpapar.

5.2.2 Odds Ratio (OR)

Odds Ratio membandingkan peluang kejadian antara dua kelompok, berguna terutama dalam studi kasus-kontrol.

Rumus :

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

dengan

$a = 5 (Beli, tanpa logo)$
$b = 44 (Tidak beli, tanpa logo)$
$c = 0 (Beli, dengan logo)$
$d = 96(Tidak beli, dengan logo)$

a <- 5
b <- 44
c <- 0
d <- 96

OR <- (a * d) / (b * c)
OR

## [1] Inf

Hasil dan Interpretasi

\[ OR = \frac{5 \times 96}{44 \times 0} = \infty \]

Sama seperti RR, nilai OR tidak terdefinisi karena tidak ada kejadian pada kelompok tidak terpapar. Ini menunjukkan bahwa peluang beli hanya terjadi di kelompok tanpa logo halal.’

Tabel Ukuran Asosiasi

Ukuran	Definisi	Desain Sampling	Interpretasi
Risk Difference (RD)	Selisih risiko antara dua kelompok. \[RD = P_{ekspos} - P_{tidak\_ekspos}\]	Kohort, eksperimen acak, eksperimen klinis	Perbedaan risiko absolut. RD = 0 berarti tidak ada perbedaan.
Relative Risk (RR)	Perbandingan risiko antara dua kelompok. \[RR = \frac{P_{ekspos}}{P_{tidak\_ekspos}}\]	Kohort, eksperimen acak, eksperimen klinis	RR > 1 artinya risiko lebih tinggi pada kelompok terpapar.
Odds Ratio (OR)	Perbandingan odds kejadian antara dua kelompok. \[OR = \frac{ad}{bc}\]	Kasus-kontrol, eksperimen klinis, observasional	OR > 1 artinya lebih besar kemungkinan kejadian pada terpapar.

Kesimpulan

Ketiga ukuran statistik, Risk Difference (RD), Relative Risk (RR), dan Odds Ratio (OR), digunakan untuk mengukur hubungan antara paparan dan kejadian. RD mengukur perbedaan risiko absolut, RR mengukur perbandingan risiko relatif, dan OR mengukur perbandingan odds kejadian. RD dan RR cocok untuk studi kohort dan eksperimen, sementara OR lebih sering digunakan dalam studi kasus-kontrol. Masing-masing memberikan informasi berbeda tentang risiko dan kejadian.

6. Inferensi Tabel Kontigensi Dua Arah

Inferensi tabel kontingensi dua arah digunakan untuk menguji hubungan antara dua variabel kategorikal. Tabel kontingensi menyajikan frekuensi kejadian dari kombinasi kategori kedua variabel tersebut. Tujuan inferensi adalah untuk menentukan apakah ada asosiasi yang signifikan antara variabel-variabel tersebut, dengan menggunakan uji statistik seperti Chi-Square

6.1 Estimasi

Estimasi merupakan proses penarikan kesimpulan tentang parameter populasi berdasarkan data sampel.

6.1.1 Estimasi Titik

Estimasi titik adalah satu nilai angka yang digunakan sebagai perkiraan terbaik dari parameter populasi

Rumus :

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

Di mana:

$\hat{p}$ adalah proporsi sampel,
$x$ adalah jumlah keberhasilan,

$n$ adalah ukuran sampel.

6.1.2 Estimasi Interval

Estimasi interval memberikan rentang nilai yang diperkirakan memuat nilai parameter populasi dengan tingkat keyakinan tertentu, biasanya 95%

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

dimana :

$Z_{\alpha/2}$ = nilai dari distribusi normal standar untuk tingkat kepercayaan tertentu.
$\hat{p}$ = Estimasi Titik Proporsi.
$n$ = Ukuran sampel.

6.2 Uji Hipotesis

6.2.1 Uji Proporsi

Tujuan uji proporsi adalah untuk menguji apakah proporsi suatu kejadian dalam satu atau lebih kelompok berbeda secara signifikan dari suatu nilai tertentu atau antar kelompok.

	Kejadian (+)	Tidak Kejadian (−)	Total
Kelompok 1	\[n_{11}\]	\[n_{12}\]	\[n_{1.}\]
Kelompok 2	\[n_{21}\]	\[n_{22}\]	\[n_2.\]
Total	\[n_{.1}\]	\[n_{.2}\]	\[n\]

Tujuan dan Hipotesis

Uji proporsi dua kelompok bertujuan untuk mengetahui apakah terdapat perbedaan proporsi kejadian antara dua kelompok. Uji ini dilakukan dengan menggunakan uji statistik z dua proporsi, dengan formulasi hipotesis sebagai berikut:

Hipotesis Nol (H0)

Proporsi kedua kelompok sama, yaitu \[p_1 = p_2\]
Hipotesis Alternatif (H1):
Proporsi kedua kelompok berbeda, yaitu \[p_1 \ne p_2\]

Estimasi Proporsi

Proporsi kejadian di kelompok 1: \[\hat{p}1 = \frac{n{11}}{n_{1.}}\]
Proporsi kejadian di kelompok 2: \[\hat{p}2 = \frac{n{21}}{n_{2.}}\]

Proporsi Gabungan

Jika diasumsikan bahwa H0 benar, maka digunakan proporsi gabungan sebagai:

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

Statistik Uji

Nilai statistik z dihitung dengan

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

Keputusan Uji

Statistik Z mengikuti distribusi normal standar \[N(0,1)\].
Bandingkan nilai |Z|

dengan nilai kritis \[Z_{\alpha/2}\] contoh: 1.96 untuk \[\alpha = 0.05\].
Jika |Z| lebih besar dari nilai kritis, maka tolak H0 → artinya terdapat perbedaan proporsi yang signifikan antara dua kelompok.

Langkah 1: Hitung Proporsi Sampel

Kelompok 1 (tidak ada logo halal):
\[n_{11} = 5\], \[n_{1.} = 49\]
\[\Rightarrow \hat{p}_1 = \frac{5}{49} \approx 0.1020\]
Kelompok 2 (ada logo halal):
\[n_{21} = 0\], \[n_{2.} = 96\]
\[\Rightarrow \hat{p}_2 = \frac{0}{96} = 0\]

Langkah 2: Hitung Proporsi Gabungan

Gabungan: \[n_{11} + n_{21} = 5 + 0 = 5\]
\[n_{1.} + n_{2.} = 49 + 96 = 145\]

\[\Rightarrow \hat{p} = \frac{5}{145} \approx 0.0345\]

Langkah 3: Hitung Statistik Uji Z

\[ Z = \frac{0.1020 - 0.0000} {\sqrt{0.0345(1 - 0.0345)\left(\frac{1}{49} + \frac{1}{96}\right)}} \]

\[ Z = \frac{0.1020} {\sqrt{0.0345 \times 0.9655 \times (0.0204 + 0.0104)}} \]

\[ Z = \frac{0.1020} {\sqrt{0.0333 \times 0.0308}} = \frac{0.1020}{\sqrt{0.001026}} = \frac{0.1020}{0.0320} \approx 3.18 \]

Interpretasi :

Jika \[Z \approx 3.18\] dan nilai kritis Z untuk \[\alpha = 0.05\] adalah 1.96, maka:

Karena |Z| > 1.96, tolak hipotesis nol sehingga terdapat perbedaan signifikan proporsi pembelian antara kelompok yang tidak ada logo halal dan yang ada.

Perhitungan dengan software R :

p1 <- 5 / 49
p2 <- 0 / 96
p1

## [1] 0.1020408

p2

## [1] 0

p_pool <- 5 / 145
p_pool

## [1] 0.03448276

n1 <- 49
n2 <- 96
p_diff <- p1 - p2
se <- sqrt(p_pool * (1 - p_pool) * (1/n1 + 1/n2))
z <- p_diff / se
z

## [1] 3.185243

6.2.2 Uji Asosiasi

Tujuan uji asosiasi dalam tabel kontingensi 2×2 adalah untuk mengetahui apakah terdapat hubungan yang signifikan antara dua variabel kategorik.

1. Risk Difference (RD)

Mengukur perbedaan absolut dari risiko kejadian antara dua kelompok.

Rumus :

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

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

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

2. Relative Risk (RR)

Mengukur seberapa besar kemungkinan kejadian pada satu kelompok dibandingkan kelompok lain.

Rumus :

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

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

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

3. Odds Ratio (OR)

Mengukur perbandingan odds antara dua kelompok.

Rumus :

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

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

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

Contoh Perhitungan Manual

Untuk Risk Difference (RD)

Rumus RD:

\[ RD = P_1 - P_2 \]

\[ = \frac{5}{49} - \frac{0}{96} = 0.1020 - 0 = 0.1020 \]

Standard Error RD:

\[ SE_{RD} = \sqrt{ \frac{p_1 (1 - p_1)}{n_1} + \frac{p_2 (1 - p_2)}{n_2} } \]

\[ = \sqrt{ \frac{0.1020 (1 - 0.1020)}{49} + \frac{0 (1 - 0)}{96} } \]

\[ = \sqrt{ \frac{0.0916}{49} } = \sqrt{0.00187} = 0.0432 \]

Z statistik:

\[ Z = \frac{RD}{SE_{RD}} = \frac{0.1020}{0.0432} \approx 2.36 \]

Untuk Relative Risk (RR)

Rumus RR:

\[ RR = \frac{p_1}{p_2} = \frac{0.1020}{0} \]

→ Tidak terdefinisi karena $p_2 = 0$ (pembilang ada, penyebut nol).

Tapi jika kita gunakan continuity correction dengan 0.5:

$p_1 = \frac{5 + 0.5}{49 + 1} = \frac{5.5}{50} = 0.11$
$p_2 = \frac{0 + 0.5}{96 + 1} = \frac{0.5}{97} \approx 0.0052$

\[ RR = \frac{0.11}{0.0052} \approx 21.15 \]

SE untuk log(RR):

\[ SE_{\log(RR)} = \sqrt{ \frac{1}{a} - \frac{1}{a + b} + \frac{1}{c} - \frac{1}{c + d} } \]

(a = 5, b = 44, c = 0, d = 96 → c = 0, tak terdefinisi)

→ Maka pakai continuity correction: a = 5.5, b = 44.5, c = 0.5, d = 96.5

\[ SE_{\log(RR)} = \sqrt{ \frac{1}{5.5} + 0.5 } = \sqrt{0.182 + 2} = \sqrt{2.182} \approx 1.477 \]

Z statistik:

\[ Z = \frac{\log(RR)}{SE_{\log(RR)}} = \frac{\log(21.15)}{1.477} = \frac{3.05}{1.477} \approx 2.06 \]

Untuk Odds Ratio (OR) :

OR (pakai continuity correction):

\[ OR = \frac{(a \cdot d)}{(b \cdot c)} = \frac{(5.5 \cdot 96.5)}{(44.5 \cdot 0.5)} \approx \frac{530.75}{22.25} \approx 23.86 \]

SE untuk log(OR):

\[ SE_{\log(OR)} = \sqrt{ \frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d} } \]

\[ = \sqrt{ \frac{1}{5.5} + \frac{1}{44.5} + \frac{1}{0.5} + \frac{1}{96.5} } \approx \sqrt{0.182 + 0.022 + 2 + 0.010} = \sqrt{2.214} \approx 1.488 \]

Z statistik:

\[ Z = \frac{\log(OR)}{SE_{\log(OR)}} = \frac{\log(23.86)}{1.488} = \frac{3.17}{1.488} \approx 2.13 \]

# Data
a <- 5     # Beli, Tidak ada logo halal
b <- 44    # Tidak beli, Tidak ada logo halal
c <- 0     # Beli, Ada logo halal
d <- 96    # Tidak beli, Ada logo halal

# Tambahkan continuity correction untuk RR dan OR (menghindari pembagi nol)
a_cc <- a + 0.5
b_cc <- b + 0.5
c_cc <- c + 0.5
d_cc <- d + 0.5

n1 <- a + b   # Tidak ada logo halal
n2 <- c + d   # Ada logo halal

# ---------------------------------------
# 1. Risk Difference (RD)
p1 <- a / n1
p2 <- c / n2
RD <- p1 - p2

# Standard error RD
se_RD <- sqrt((p1 * (1 - p1)) / n1 + (p2 * (1 - p2)) / n2)

# Z statistic for RD
z_RD <- RD / se_RD

# ---------------------------------------
# 2. Relative Risk (RR)
RR <- (a_cc / n1) / (c_cc / n2)

# Standard error for log(RR)
se_logRR <- sqrt((1 / a_cc) - (1 / n1) + (1 / c_cc) - (1 / n2))
z_RR <- log(RR) / se_logRR

# ---------------------------------------
# 3. Odds Ratio (OR)
OR <- (a_cc * d_cc) / (b_cc * c_cc)

# Standard error for log(OR)
se_logOR <- sqrt(1 / a_cc + 1 / b_cc + 1 / c_cc + 1 / d_cc)
z_OR <- log(OR) / se_logOR

# ---------------------------------------
# Output
cat("Risk Difference (RD):", round(RD, 4), "\n")

## Risk Difference (RD): 0.102

cat("Z statistic (RD):", round(z_RD, 2), "\n\n")

## Z statistic (RD): 2.36

cat("Relative Risk (RR):", round(RR, 2), "\n")

## Relative Risk (RR): 21.55

cat("Z statistic (RR):", round(z_RR, 2), "\n\n")

## Z statistic (RR): 2.09

cat("Odds Ratio (OR):", round(OR, 2), "\n")

## Odds Ratio (OR): 23.85

cat("Z statistic (OR):", round(z_OR, 2), "\n")

## Z statistic (OR): 2.13

6.2.3 Uji Independensi

Tujuan uji independensi adalah untuk mengetahui apakah terdapat hubungan atau ketergantungan antara dua variabel kategori dalam suatu populasi.

6.2.3.1 Uji Chi-Square

Tujuan uji Chi-Square adalah untuk menguji apakah terdapat hubungan yang signifikan antara dua variabel kategori dalam suatu populasi.

Rumus :

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

Keterangan:

$O_{ij}$: Frekuensi observasi pada sel ke-$i, j$
$E_{ij}$: Frekuensi ekspektasi (harapan) pada sel ke-$i, j$
$\sum$: Penjumlahan atas semua sel dalam tabel

Perhitungan manual :

1. Hitung Frekuensi Harapan (E) untuk masing-masing sel:

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

Contoh:

$E_{11} = \frac{5 \times 49}{145} = 1.69$
$E_{12} = \frac{5 \times 96}{145} = 3.31$
$E_{21} = \frac{140 \times 49}{145} = 47.31$
$E_{22} = \frac{140 \times 96}{145} = 92.69$

2. Hitung nilai Chi-Square:

\[ \chi^2 = \frac{(5 - 1.69)^2}{1.69} + \frac{(0 - 3.31)^2}{3.31} + \frac{(44 - 47.31)^2}{47.31} + \frac{(96 - 92.69)^2}{92.69} \]

Langkah-langkah:

\[ \chi^2 = \frac{3.31^2}{1.69} + \frac{3.31^2}{3.31} + \frac{(-3.31)^2}{47.31} + \frac{3.31^2}{92.69} \]

\[ \chi^2 \approx \frac{10.96}{1.69} + \frac{10.96}{3.31} + \frac{10.96}{47.31} + \frac{10.96}{92.69} \]

\[ \chi^2 \approx 6.49 + 3.31 + 0.23 + 0.12 = \boxed{10.15} \]

Perhitungan dengan software R :

# Buat matriks kontingensi
tabel <- matrix(c(5, 0, 44, 96),
                nrow = 2,
                byrow = TRUE)

# Tambahkan nama baris dan kolom
dimnames(tabel) <- list(
  Keputusan = c("Beli", "Tidak Beli"),
  LogoHalal = c("Tidak Ada", "Ada")
)

# Lihat tabel
tabel

##             LogoHalal
## Keputusan    Tidak Ada Ada
##   Beli               5   0
##   Tidak Beli        44  96

# Uji chi-square
chisq.test(tabel, correct = FALSE)

## Warning in chisq.test(tabel, correct = FALSE): Chi-squared approximation may be
## incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tabel
## X-squared = 10.146, df = 1, p-value = 0.001446

Interpretasi :

Berdasarkan hasil uji Chi-Square dengan nilai statistik sebesar 10.15 dan derajat kebebasan 1, nilai tersebut melebihi nilai kritis pada taraf signifikansi 5% (yakni sekitar 3.84). Hal ini menunjukkan bahwa terdapat hubungan yang signifikan antara keberadaan logo halal pada produk dan keputusan konsumen untuk membeli. Dengan kata lain, keputusan membeli konsumen dipengaruhi oleh ada tidaknya logo halal pada produk.

6.2.3.2 Uji Likelihood Ratio

\[ G^2 = 2 \sum O_{ij} \cdot \ln\left( \frac{O_{ij}}{E_{ij}} \right) \]

Keterangan:

$O_{ij}$: Frekuensi observasi pada sel ke-$i, j$
$E_{ij}$: Frekuensi ekspektasi (expected) pada sel ke-$i, j$
$\ln$: Logaritma natural (basis $e$)

Perhitungan manual :

1. Hitung nilai harapan (E)

\[ E_{11} = \frac{5 \times 49}{145} = 1.69, \quad E_{12} = \frac{5 \times 96}{145} = 3.31 \]

\[ E_{21} = \frac{140 \times 49}{145} = 47.31, \quad E_{22} = \frac{140 \times 96}{145} = 92.69 \]

2. Substitusi ke rumus

\[ G^2 = 2 \times \left[ 5 \cdot \ln\left( \frac{5}{1.69} \right) + 0 \cdot \ln\left( \frac{0}{3.31} \right) + 44 \cdot \ln\left( \frac{44}{47.31} \right) + 96 \cdot \ln\left( \frac{96}{92.69} \right) \right] \]

\[ G^2 = 2 \times [5 \cdot \ln(2.96) + 0 + 44 \cdot \ln(0.93) + 96 \cdot \ln(1.04)] \]

\[ G^2 = 2 \times [5 \cdot 1.08 + 44 \cdot (-0.07) + 96 \cdot 0.04] \]

\[ G^2 = 2 \times [5.4 - 3.08 + 3.84] = 2 \times 6.16 = \mathbf{12.32} \]

Perhitungan dengan software R :

# Load package
library(MASS)

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

# Buat tabel kontingensi
tabel <- matrix(c(5, 0, 44, 96),
                nrow = 2,
                byrow = TRUE)

# Tambahkan nama dimensi
dimnames(tabel) <- list(
  Keputusan = c("Beli", "Tidak Beli"),
  LogoHalal = c("Tidak Ada", "Ada")
)

# Lakukan uji log-likelihood ratio
model_loglm <- loglm(~ Keputusan + LogoHalal, data = tabel)

# Lihat hasil
summary(model_loglm)

## Formula:
## ~Keputusan + LogoHalal
## attr(,"variables")
## list(Keputusan, LogoHalal)
## attr(,"factors")
##           Keputusan LogoHalal
## Keputusan         1         0
## LogoHalal         0         1
## attr(,"term.labels")
## [1] "Keputusan" "LogoHalal"
## 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 11.20321  1 0.0008165615
## Pearson          10.14577  1 0.0014463253

Interpretasi :

Karena p-value < 0.05, maka kita menolak H₀ yang menyatakan bahwa variabel “Keputusan Membeli” dan “Logo Halal” saling bebas atau tidak berasosiasi.

6.2.3.3 Uji Exact Fisher

Uji Exact Fisher digunakan untuk menganalisis tabel kontingensi 2x2, terutama bila ukuran sampel kecil atau ada nilai frekuensi yang sangat kecil (bahkan nol). Uji ini menghitung probabilitas mendapatkan distribusi frekuensi tertentu (atau yang lebih ekstrem) dengan asumsi tidak ada asosiasi antara dua variabel.

\[ P = \frac{(a + b)! \cdot (c + d)! \cdot (a + c)! \cdot (b + d)!}{a! \cdot b! \cdot c! \cdot d! \cdot n!} \]

Keterangan:

$a, b, c, d$: frekuensi pada sel
$!$: faktorial
$n$: total observasi

Perhitungan manual :

\[ P = \frac{(5 + 0)! \cdot (44 + 96)! \cdot (5 + 44)! \cdot (0 + 96)!}{5! \cdot 0! \cdot 44! \cdot 96! \cdot 145!} \]

Perhitungan menggunakan software R :

# Membuat matriks kontingensi
data <- matrix(c(5, 0, 44, 96), nrow = 2, byrow = TRUE)

# Menambahkan nama baris dan kolom
dimnames(data) <- list(
  Keputusan = c("Beli", "Tidak Beli"),
  LogoHalal = c("Tidak Ada", "Ada")
)

# Tampilkan tabel
print(data)

##             LogoHalal
## Keputusan    Tidak Ada Ada
##   Beli               5   0
##   Tidak Beli        44  96

# Lakukan uji Fisher
fisher.test(data)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  data
## p-value = 0.003828
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.890899      Inf
## sample estimates:
## odds ratio 
##        Inf

Interpretasi :

Karena p-value < 0.05, maka kita menolak H₀ yang menyatakan bahwa ada asosiasi signifikan antara keberadaan logo halal dan keputusan beli

6.3 Analisis Residual dalam Tabel Kontigensi

Analisis residual digunakan untuk mengetahui sel mana dalam tabel kontingensi yang memberikan kontribusi paling besar terhadap hasil uji chi-square, yaitu deviasi antara frekuensi yang diamati (observed) dan yang diharapkan (expected) berdasarkan hipotesis nol. Analisis ini sangat berguna untuk mengidentifikasi pola atau hubungan yang mungkin tersembunyi di dalam data.

Pearson Residual

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

Keterangan:

$R_{ij}$: nilai residual Pearson pada sel ke-$i,j$
$O_{ij}$: nilai observasi (frekuensi aktual)
$E_{ij}$: nilai harapan (frekuensi yang diharapkan berdasarkan hipotesis nol)

Perhitungan Manual:

1. Menghitung Frekuensi Harapan

$E_{11} = \frac{5 \cdot 49}{145} \approx 1.69$
$E_{12} = \frac{5 \cdot 96}{145} \approx 3.31$
$E_{21} = \frac{140 \cdot 49}{145} \approx 47.31$
$E_{22} = \frac{140 \cdot 96}{145} \approx 92.69$

2. Menghitung Residual Pearson untuk setiap sel

Untuk sel Beli & Tidak Ada Logo Halal:
\[ R = \frac{5 - 1.69}{\sqrt{1.69}} \approx 2.55 \]
Untuk sel Beli & Ada Logo Halal:
\[ R = \frac{0 - 3.31}{\sqrt{3.31}} \approx -1.82 \]
Untuk sel Tidak Beli & Tidak Ada Logo Halal:
\[ R = \frac{44 - 47.31}{\sqrt{47.31}} \approx -0.48 \]
Untuk sel Tidak Beli & Ada Logo Halal:
\[ R = \frac{96 - 92.69}{\sqrt{92.69}} \approx 0.34 \]

Perhitungan menggunakan software R :

# Membuat matriks kontingensi
data <- matrix(c(5, 0, 44, 96), nrow = 2, byrow = TRUE)
dimnames(data) <- list(
  Keputusan = c("Beli", "Tidak Beli"),
  LogoHalal = c("Tidak Ada", "Ada")
)

# Melakukan uji chi-square (tanpa continuity correction)
uji_chi <- chisq.test(data, correct = FALSE)

## Warning in chisq.test(data, correct = FALSE): Chi-squared approximation may be
## incorrect

# Menampilkan residual Pearson
uji_chi$residuals

##             LogoHalal
## Keputusan     Tidak Ada       Ada
##   Beli        2.5466789 -1.819435
##   Tidak Beli -0.4812771  0.343841

Kesimpulan

Tidak ada Outlier karena tidak ada nilai |r| > 3
Sel “Beli & Tidak Ada Logo Halal” merupakan satu-satunya sel yang secara signifikan menyimpang dari nilai harapan, dengan residual +2.55.

7. Tabel Kontigensi Tiga Arah

Tabel kontingensi tiga arah adalah bentuk perluasan dari tabel kontingensi dua arah yang digunakan untuk menyajikan dan menganalisis hubungan tiga variabel kategorik secara simultan.

Jika tabel dua arah menyajikan hubungan antara dua variabel kategorik (misalnya, jenis kelamin dan status merokok), maka tabel tiga arah menambahkan dimensi ketiga yang dapat mempengaruhi atau mengubah hubungan antara dua variabel awal.

Tujuan utama adalah untuk:

Mengetahui apakah hubungan antara dua variabel (misalnya A dan B) konsisten di seluruh kategori variabel ketiga (C).
Mengidentifikasi adanya interaksi atau efek pengendalian (control effect) oleh variabel ketiga.

7.1 Tabel Parsial dan Marginal

Tabel Parsial adalah tabel dua arah yang disusun dengan mengendalikan atau memfiksasi salah satu variabel dari tabel kontingensi tiga arah.

Tabel Marginal adalah tabel dua arah yang dihasilkan dengan mengabaikan atau menggabungkan (menjumlahkan) dimensi dari salah satu variabel dalam tabel tiga arah.

Contoh data :

Jenis Kelamin	Pendidikan	Setuju	Tidak Setuju
Pria	≤ 8	72	47
Pria	9 – 12	110	196
Pria	≥ 13	44	179
Wanita	≤ 8	158	85
Wanita	9 – 12	173	283
Wanita	≥ 13	28	187

Tabel Frekuensi Parsial

# Membuat array 3 dimensi 
data <- array(
  c(
    72, 47,    # Pria, <=8
    110, 196,  # Pria, 9–12
    44, 179,   # Pria, >=13
    158, 85,   # Wanita, <=8
    173, 283,  # Wanita, 9–12
    28, 187    # Wanita, >=13
  ),
  dim = c(2, 3, 2),
  dimnames = list(
    Pendapat = c("Setuju", "Tidak setuju"),
    Pendidikan = c("<=8", "9-12", ">=13"),
    JenisKelamin = c("Pria", "Wanita")
  )
)

  # Tabel Parsial untuk Pria
data[, , "Pria"]

##               Pendidikan
## Pendapat       <=8 9-12 >=13
##   Setuju        72  110   44
##   Tidak setuju  47  196  179

# Tabel Parsial untuk Wanita
data[, , "Wanita"]

##               Pendidikan
## Pendapat       <=8 9-12 >=13
##   Setuju       158  173   28
##   Tidak setuju  85  283  187

Tabel Frekuensi Marginal

# Tabel data 3 dimensi sudah dibuat sebelumnya sebagai `data`

# Total keseluruhan
total <- sum(data)

# Marginal menurut Pendapat & Pendidikan (menjumlahkan per Jenis Kelamin)
margin_pendapat_pendidikan <- apply(data, c(1, 2), sum)
print("Marginal: Pendapat x Pendidikan")

## [1] "Marginal: Pendapat x Pendidikan"

print(margin_pendapat_pendidikan)

##               Pendidikan
## Pendapat       <=8 9-12 >=13
##   Setuju       230  283   72
##   Tidak setuju 132  479  366

# Marginal menurut Pendapat & Jenis Kelamin (menjumlahkan per Pendidikan)
margin_pendapat_gender <- apply(data, c(1, 3), sum)
print("Marginal: Pendapat x Jenis Kelamin")

## [1] "Marginal: Pendapat x Jenis Kelamin"

print(margin_pendapat_gender)

##               JenisKelamin
## Pendapat       Pria Wanita
##   Setuju        226    359
##   Tidak setuju  422    555

# Marginal menurut Pendidikan & Jenis Kelamin (menjumlahkan per Pendapat)
margin_pendidikan_gender <- apply(data, c(2, 3), sum)
print("Marginal: Pendidikan x Jenis Kelamin")

## [1] "Marginal: Pendidikan x Jenis Kelamin"

print(margin_pendidikan_gender)

##           JenisKelamin
## Pendidikan Pria Wanita
##       <=8   119    243
##       9-12  306    456
##       >=13  223    215

# Marginal tunggal (jumlah tiap dimensi)
margin_pendapat <- apply(data, 1, sum)
margin_pendidikan <- apply(data, 2, sum)
margin_gender <- apply(data, 3, sum)

print("Marginal: Pendapat")

## [1] "Marginal: Pendapat"

print(margin_pendapat)

##       Setuju Tidak setuju 
##          585          977

print("Marginal: Pendidikan")

## [1] "Marginal: Pendidikan"

print(margin_pendidikan)

##  <=8 9-12 >=13 
##  362  762  438

print("Marginal: Jenis Kelamin")

## [1] "Marginal: Jenis Kelamin"

print(margin_gender)

##   Pria Wanita 
##    648    914

7.2 Distribusi Peluang

1. Peluang Bersama (Joint Probability)

\[ P(A, B, C) = \frac{f(A, B, C)}{N} \]

Keterangan:

f(A, B, C): frekuensi kasus gabungan A, B, dan C
N: total seluruh frekuensi

# Peluang bersama: setiap sel dibagi total keseluruhan
joint_prob <- round(data / sum(data), 4)
joint_prob

## , , JenisKelamin = Pria
## 
##               Pendidikan
## Pendapat          <=8   9-12   >=13
##   Setuju       0.0461 0.0704 0.0282
##   Tidak setuju 0.0301 0.1255 0.1146
## 
## , , JenisKelamin = Wanita
## 
##               Pendidikan
## Pendapat          <=8   9-12   >=13
##   Setuju       0.1012 0.1108 0.0179
##   Tidak setuju 0.0544 0.1812 0.1197

2. Peluang Marginal

\[ P(A) = \frac{\sum_{b,c} f(A, B, C)}{N} \]

Keterangan:

Menjumlahkan seluruh nilai pada dimensi lainnya
Digunakan untuk melihat distribusi masing-masing kategori secara umum

# Peluang marginal Pendapat
marg_pendapat <- round(apply(data, 1, sum) / sum(data), 4)
print("Peluang marginal Pendapat:")

## [1] "Peluang marginal Pendapat:"

print(marg_pendapat)

##       Setuju Tidak setuju 
##       0.3745       0.6255

# Peluang marginal Pendidikan
marg_pendidikan <- round(apply(data, 2, sum) / sum(data), 4)
print("Peluang marginal Pendidikan:")

## [1] "Peluang marginal Pendidikan:"

print(marg_pendidikan)

##    <=8   9-12   >=13 
## 0.2318 0.4878 0.2804

# Peluang marginal Jenis Kelamin
marg_gender <- round(apply(data, 3, sum) / sum(data), 4)
print("Peluang marginal Jenis Kelamin:")

## [1] "Peluang marginal Jenis Kelamin:"

print(marg_gender)

##   Pria Wanita 
## 0.4149 0.5851

3. Peluang Bersyarat

\[ P(A \mid B) = \frac{P(A, B)}{P(B)} = \frac{f(A, B)}{\sum_{a} f(A, B)} \]

Keterangan:

Menunjukkan peluang A terjadi dengan syarat B sudah terjadi

# Peluang bersyarat: Pendapat | Jenis Kelamin = Pria
pria_total <- apply(data[, , "Pria"], 1, sum)
prob_setuju_given_pria <- data["Setuju", , "Pria"] / sum(data[, , "Pria"])
print("Peluang Setuju untuk tiap tingkat pendidikan (diberi Pria):")

## [1] "Peluang Setuju untuk tiap tingkat pendidikan (diberi Pria):"

print(round(prob_setuju_given_pria, 4))

##    <=8   9-12   >=13 
## 0.1111 0.1698 0.0679

# Peluang bersyarat: Pendapat | Pendidikan = "9-12"
pend_912_total <- sum(data[, "9-12", ])
prob_setuju_given_pend912 <- data["Setuju", "9-12", ] / pend_912_total
print("Peluang Setuju (diberi Pendidikan 9-12):")

## [1] "Peluang Setuju (diberi Pendidikan 9-12):"

print(round(prob_setuju_given_pend912, 4))

##   Pria Wanita 
## 0.1444 0.2270

# Peluang bersyarat: Pendidikan | Pendapat = Setuju
pendapat_setuju_total <- sum(data["Setuju", , ])
prob_pendidikan_given_setuju <- data["Setuju", , ] / pendapat_setuju_total
print("Peluang Pendidikan (diberi Setuju):")

## [1] "Peluang Pendidikan (diberi Setuju):"

print(round(prob_pendidikan_given_setuju, 4))

##           JenisKelamin
## Pendidikan   Pria Wanita
##       <=8  0.1231 0.2701
##       9-12 0.1880 0.2957
##       >=13 0.0752 0.0479

7.3 Tabel Peluang Bersyarat

# Buat array 3 dimensi (jika belum dibuat sebelumnya)
data <- array(
  c(
    72, 47,    # Pria, ≤8
    110, 196,  # Pria, 9–12
    44, 179,   # Pria, ≥13
    158, 85,   # Wanita, ≤8
    173, 283,  # Wanita, 9–12
    28, 187    # Wanita, ≥13
  ),
  dim = c(2, 3, 2),
  dimnames = list(
    Pendapat = c("Setuju", "Tidak setuju"),
    Pendidikan = c("<=8", "9-12", ">=13"),
    JenisKelamin = c("Pria", "Wanita")
  )
)

# Hitung peluang bersyarat berdasarkan Jenis Kelamin
# Ini akan menghasilkan array 2D untuk setiap Jenis Kelamin
prob_given_gender <- prop.table(data, margin = 3)
round(prob_given_gender, 4)  # Dibulatkan 4 angka di belakang koma

## , , JenisKelamin = Pria
## 
##               Pendidikan
## Pendapat          <=8   9-12   >=13
##   Setuju       0.1111 0.1698 0.0679
##   Tidak setuju 0.0725 0.3025 0.2762
## 
## , , JenisKelamin = Wanita
## 
##               Pendidikan
## Pendapat          <=8   9-12   >=13
##   Setuju       0.1729 0.1893 0.0306
##   Tidak setuju 0.0930 0.3096 0.2046

7.4 Ukuran Asosiasi

Risk Difference

Selisih risiko antara kelompok terpapar dan tidak terpapar.

Relative Risk

Perbandingan risiko kejadian antara kelompok terpapar dan tidak terpapar.

Odds Ratio

Perbandingan odds (peluang dibanding bukan peluang) antara dua kelompok.

# Data array 3D
data <- array(
  c(
    72, 47,    # Pria, ≤8
    110, 196,  # Pria, 9–12
    44, 179,   # Pria, ≥13
    158, 85,   # Wanita, ≤8
    173, 283,  # Wanita, 9–12
    28, 187    # Wanita, ≥13
  ),
  dim = c(2, 3, 2),
  dimnames = list(
    Pendapat = c("Setuju", "Tidak setuju"),
    Pendidikan = c("<=8", "9-12", ">=13"),
    JenisKelamin = c("Pria", "Wanita")
  )
)

# Menjumlahkan semua tingkat pendidikan -> hasil tabel 2x2
tabel_bersama <- apply(data, c(1, 3), sum)
tabel_bersama

##               JenisKelamin
## Pendapat       Pria Wanita
##   Setuju        226    359
##   Tidak setuju  422    555

# Ekstrak nilai
a <- tabel_bersama["Setuju", "Pria"]
b <- tabel_bersama["Tidak setuju", "Pria"]
c <- tabel_bersama["Setuju", "Wanita"]
d <- tabel_bersama["Tidak setuju", "Wanita"]

# Hitung risiko
risk_pria <- a / (a + b)
risk_wanita <- c / (c + d)

# Risk Difference
RD <- risk_wanita - risk_pria

# Relative Risk
RR <- risk_wanita / risk_pria

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

# Cetak hasil
cat("Risk Difference (RD):", round(RD, 4), "\n")

## Risk Difference (RD): 0.044

cat("Relative Risk (RR):", round(RR, 4), "\n")

## Relative Risk (RR): 1.1262

cat("Odds Ratio (OR):", round(OR, 4), "\n")

## Odds Ratio (OR): 1.2078

Interpretasi :

RD: Wanita memiliki proporsi persetujuan 4.3% lebih tinggi dibanding pria secara keseluruhan.
RR: Wanita 1.1 kali lebih mungkin menyatakan setuju dibanding pria.
OR: Peluang wanita menyatakan setuju sekitar 1.11 kali peluang pria secara keseluruhan.

7.4.1 Tabel Kontigensi Parsial

# Data array 3D
data <- array(
  c(
    72, 47,    # Pria, ≤8
    110, 196,  # Pria, 9–12
    44, 179,   # Pria, ≥13
    158, 85,   # Wanita, ≤8
    173, 283,  # Wanita, 9–12
    28, 187    # Wanita, ≥13
  ),
  dim = c(2, 3, 2),
  dimnames = list(
    Pendapat = c("Setuju", "Tidak setuju"),
    Pendidikan = c("<=8", "9-12", ">=13"),
    JenisKelamin = c("Pria", "Wanita")
  )
)

# Ambil data untuk Pendidikan = ">=13"
# Baris = Pendapat, Kolom = Jenis Kelamin
tabel_2x2 <- data[, ">=13", ]
tabel_2x2

##               JenisKelamin
## Pendapat       Pria Wanita
##   Setuju         44     28
##   Tidak setuju  179    187

# Ambil nilai dari tabel
a <- tabel_2x2["Setuju", "Pria"]
b <- tabel_2x2["Tidak setuju", "Pria"]
c <- tabel_2x2["Setuju", "Wanita"]
d <- tabel_2x2["Tidak setuju", "Wanita"]

# Hitung risiko
risk_pria <- a / (a + b)
risk_wanita <- c / (c + d)

# RD
RD <- risk_wanita - risk_pria

# RR
RR <- risk_wanita / risk_pria

# OR
OR <- (c / d) / (a / b)

# Tampilkan hasil
cat("Risk Difference (RD):", round(RD, 4), "\n")

## Risk Difference (RD): -0.0671

cat("Relative Risk (RR):", round(RR, 4), "\n")

## Relative Risk (RR): 0.66

cat("Odds Ratio (OR):", round(OR, 4), "\n")

## Odds Ratio (OR): 0.6091

Interpretasi :

RD: Proporsi wanita yang setuju 13.4% lebih rendah dibanding pria pada tingkat pendidikan ≥13.

RR: Wanita 0.68 kali lebih mungkin menyatakan setuju dibanding pria.

OR: Peluang wanita menyatakan setuju 0.65 kali dibanding pria pada tingkat pendidikan ≥13.

7.5 Conditional Independence

Conditional Independence (Independensi Bersyarat) terjadi ketika dua variabel (misalnya A dan B) tidak saling berhubungan, setelah kita mengontrol atau menahan variabel ketiga (C) tetap konstan.

Tujuan Conditional Independence adalah apakah hubungan antara dua variabel (misalnya Pendapat dan Jenis Kelamin) masih ada atau hilang setelah dikontrol terhadap variabel ketiga (Pendidikan)

# Data array 3D
data <- array(
  c(
    72, 47,    # Pria, ≤8
    110, 196,  # Pria, 9–12
    44, 179,   # Pria, ≥13
    158, 85,   # Wanita, ≤8
    173, 283,  # Wanita, 9–12
    28, 187    # Wanita, ≥13
  ),
  dim = c(2, 3, 2),
  dimnames = list(
    Pendapat = c("Setuju", "Tidak setuju"),
    Pendidikan = c("<=8", "9-12", ">=13"),
    JenisKelamin = c("Pria", "Wanita")
  )
)

library(MASS)

# Model tanpa interaksi antara Pendapat dan Jenis Kelamin
model_cond_indep <- loglm(~ Pendapat + Pendidikan + JenisKelamin + Pendapat:Pendidikan + JenisKelamin:Pendidikan, data)

summary(model_cond_indep)

## Formula:
## ~Pendapat + Pendidikan + JenisKelamin + Pendapat:Pendidikan + 
##     JenisKelamin:Pendidikan
## attr(,"variables")
## list(Pendapat, Pendidikan, JenisKelamin)
## attr(,"factors")
##              Pendapat Pendidikan JenisKelamin Pendapat:Pendidikan
## Pendapat            1          0            0                   1
## Pendidikan          0          1            0                   1
## JenisKelamin        0          0            1                   0
##              Pendidikan:JenisKelamin
## Pendapat                           0
## Pendidikan                         1
## JenisKelamin                       1
## attr(,"term.labels")
## [1] "Pendapat"                "Pendidikan"             
## [3] "JenisKelamin"            "Pendapat:Pendidikan"    
## [5] "Pendidikan:JenisKelamin"
## attr(,"order")
## [1] 1 1 1 2 2
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                       X^2 df  P(> X^2)
## Likelihood Ratio 4.625043  3 0.2014043
## Pearson          4.599688  3 0.2035688

Interpretasi :

Karena p-value > 0.05, terima H0 sehingga tidak tidak ada bukti signifikan bahwa Pendapat dan Jenis Kelamin saling terkait setelah dikontrol oleh Pendidikan.

7.6 Inferensi Tabel Kontigensi Tiga Arah

Inferensi dalam konteks tabel kontingensi 3 arah adalah proses untuk menyimpulkan hubungan antar tiga variabel kategorik, dengan menggunakan analisis statistik

7.6.1 Independensi Bersyarat dalam Tabel Kontigensi Tiga Arah

Dalam konteks analisis data kategorik, khususnya tabel kontingensi tiga arah, independensi bersyarat (conditional independence) merujuk pada suatu kondisi di mana dua variabel saling bebas (tidak memiliki asosiasi) apabila dikondisikan terhadap variabel ketiga.

Secara formal, misalkan terdapat tiga variabel kategorik:

A = Pendapat
B = Jenis Kelamin
C = Pendidikan

Maka, $A \perp B \mid C$ menyatakan bahwa A dan B bersifat independen bersyarat terhadap C, yaitu:

Dalam setiap kategori dari variabel C, tidak terdapat hubungan (asosiasi) antara A dan B.

7.6.2 Pengujian Statistik untuk Independensi Bersyarat

Uji Cochran-Maentel-Haenzel (CMH) merupakan uji yang digunakan untuk menguji asosiasi antara dua variabel kategorik (2×2), dengan mengontrol pengaruh dari variabel ketiga (stratifikasi).

\[ \text{CMH} = \frac{\left[ \sum_{k=1}^{K} (a_k - E(a_k)) \right]^2}{\sum_{k=1}^{K} Var(a_k)} \]

Dengan:

$a_k$: frekuensi sel atas kiri (exposed + outcome positif) di strata ke-$k$
$E(a_k)$: nilai harapan dari $a_k$
$Var(a_k)$: varians dari $a_k$
CMH statistic mengikuti distribusi chi-square dengan 1 derajat kebebasan (df = 1)

# Susun array: [Pendapat, Jenis Kelamin, Pendidikan]
data <- array(c(
  # Pendidikan <=8
  72, 47,    # Pria: Setuju, Tidak setuju
  158, 85,   # Wanita

  # Pendidikan 9-12
  110, 196,
  173, 283,

  # Pendidikan >=13
  44, 179,
  28, 187
), dim = c(2, 2, 3),
dimnames = list(
  Pendapat = c("Setuju", "Tidak_setuju"),
  JenisKelamin = c("Pria", "Wanita"),
  Pendidikan = c("<=8", "9-12", ">=13")
))
mantelhaen.test(data)

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

Interpretasi :

Karena p-value > 0.05, maka terima H0 sehingga tidak terdapat asosiasi signifikan antara jenis kelamin dan pendapat jika pendidikan dikontrol ( independen secara bersyarat)

7.6.3 Odds Ratio Bersama

\[ OR_{MH} = \frac{\sum_k \frac{a_k d_k}{n_k}}{\sum_k \frac{b_k c_k}{n_k}} \]

dengan:

$k$ = indeks strata (1 sampai 3, berdasarkan pendidikan)
$n_k = a_k + b_k + c_k + d_k$

Standar Error Odds Ratio Bersama

\[ SE[\log(OR_{MH})] = \sqrt{\frac{1}{\sum_k \frac{a_k d_k}{n_k}} + \frac{1}{\sum_k \frac{b_k c_k}{n_k}}} \]

Interval Kepercayaan Odds Ratio Bersama

\[ CI_{95\%} = \exp \left[ \log(OR_{MH}) \pm Z_{1 - \alpha/2} \cdot SE \right] \]

Dengan:

$Z_{1 - \alpha/2} = 1.96$ (untuk 95%)

# Susun array 2x2x3
data <- array(c(
  72, 47,   # strata 1: Pria (Setuju, Tidak setuju)
  158, 85,  # Wanita

  110, 196,
  173, 283,

  44, 179,
  28, 187
), dim = c(2, 2, 3),
dimnames = list(
  Pendapat = c("Setuju", "Tidak_setuju"),
  JenisKelamin = c("Pria", "Wanita"),
  Pendidikan = c("<=8", "9-12", ">=13")
))


library(epitools)

## 
## Attaching package: 'epitools'

## The following object is masked from 'package:vcd':
## 
##     oddsratio

# Lakukan uji dan simpan hasilnya
hasil <- mantelhaen.test(data, correct = FALSE)
hasil

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

# Ambil nilai odds ratio bersama
or_mh <- hasil$estimate

# Ambil confidence interval
ci <- hasil$conf.int

# Hitung standar error dari log(OR)
se <- (log(ci[2]) - log(ci[1])) / (2 * 1.96)

# Tampilkan
cat("Odds Ratio Bersama (OR_MH):", or_mh, "\n")

## Odds Ratio Bersama (OR_MH): 1.001165

cat("Standard Error log(OR_MH):", se, "\n")

## Standard Error log(OR_MH): 0.1143305

cat("95% Confidence Interval:", ci[1], "sampai", ci[2], "\n")

## 95% Confidence Interval: 0.800176 sampai 1.252638

Interpretasi :

OR_MH = 1.001:
Secara keseluruhan, odds untuk menyatakan setuju nyaris sama antara pria dan wanita, setelah mengontrol pengaruh pendidikan.

95% Confidence Interval (CI) = (0.800, 1.253):
Karena CI mencakup angka 1, maka:

Tidak terdapat bukti yang cukup untuk menyimpulkan bahwa terdapat perbedaan odds yang signifikan antara pria dan wanita setelah dikontrol oleh pendidikan.

Standard Error log(OR_MH) = 0.1143:
Ini menunjukkan ketidakpastian estimasi log(OR_MH). Nilai yang relatif kecil mengindikasikan estimasi cukup stabil.

7.6.4 Uji Homogenitas Odds Ratio dengan Statistik Breslow Day

Uji Breslow-Day digunakan untuk menguji homogenitas (kesamaan) odds ratio antar beberapa strata (lapisan), misalnya untuk memastikan bahwa odds ratio pada tiap tingkat pendidikan seragam.

Rumus :

\[ Q_{BD} = \sum_{h=1}^{k} \frac{(a_h - \hat{a}_h)^2}{Var(\hat{a}_h)} \]

Dengan:

h: indeks untuk strata (misalnya tingkat pendidikan)
$a_h$: frekuensi observasi sel a (misalnya, jumlah setuju pria) pada strata ke-h
$\hat{a}_h$: frekuensi harapan berdasarkan OR gabungan (MH)
$Var(\hat{a}_h)$: varians dari frekuensi harapan

Perhitungan Manual :

Gunakan rumus:

\[ \hat{a}_1 = \frac{(a + b)(a + c)}{n} \cdot \frac{OR_{MH}}{1 + \left[(a + b)(c + d) + (a + c)(b + d)\right] \cdot \frac{(OR_{MH} - 1)}{n}} \]

Atau lebih sederhana:

\[ \hat{a}_1 = \frac{m_{1+} \cdot m_{+1} \cdot OR_{MH}}{m_{++} + (m_{1+} - m_{+1}) \cdot (OR_{MH} - 1)} \]

Dengan:

$m_{1+} = a + b = 230$
$m_{+1} = a + c = 119$
$m_{++} = n = 362$

\[ \hat{a}_1 = \frac{230 \cdot 119 \cdot 1.001165}{362 + (230 - 119) \cdot (1.001165 - 1)} \approx \frac{27326.5}{362 + 111 \cdot 0.001165} \approx \frac{27326.5}{362.129} \approx 75.47 \]

\[ \text{Var}(\hat{a}_1) = \frac{R_1 \cdot C_1 \cdot R_2 \cdot C_2}{n^2 \cdot (n - 1)} \]

Dengan:

$R_1 = a + b = 230$
$R_2 = c + d = 132$
$C_1 = a + c = 119$
$C_2 = b + d = 243$
$n = 362$

\[ \text{Var}(\hat{a}_1) = \frac{230 \cdot 119 \cdot 132 \cdot 243}{362^2 \cdot 361} = \frac{880221960}{13104484} \approx 67.16 \]

\[ Q_{BD1} = \frac{(a - \hat{a}_1)^2}{\text{Var}(\hat{a}_1)} = \frac{(72 - 75.47)^2}{67.16} = \frac{12.05}{67.16} \approx 0.179 \]

\[ Q_{BD} = Q_1 + Q_2 + Q_3 = 0.179 + 0.124 + 0.848 = \boxed{1.151} \]

library(vcdExtra)

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

## Loading required package: gnm

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

## 
## Attaching package: 'vcdExtra'

## The following object is masked from 'package:epitools':
## 
##     expand.table

library(DescTools)

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

# Data array: Pendapat (2), Pendidikan (3), JenisKelamin (2)
data <- array(
  c(
    72, 47,    # Pria, ≤8
    110, 196,  # Pria, 9–12
    44, 179,   # Pria, ≥13
    158, 85,   # Wanita, ≤8
    173, 283,  # Wanita, 9–12
    28, 187    # Wanita, ≥13
  ),
  dim = c(2, 2, 3),
  dimnames = list(
    Pendapat = c("Setuju", "TidakSetuju"),
    JenisKelamin = c("Pria", "Wanita"),
    Pendidikan = c("<=8", "9-12", ">=13")
  )
)
BreslowDayTest(data)

## 
##  Breslow-Day test on Homogeneity of Odds Ratios
## 
## data:  data
## X-squared = 161.44, df = 2, p-value < 2.2e-16

Interpretasi :

Karena p-value < 0.05, maka tolak H0 sehingga OR tidak homogen antar strata

8. Generalized Linear Model (GLM)

Generalized Linear Model (GLM) adalah perluasan dari model regresi linear klasik yang memungkinkan kita untuk menganalisis data yang tidak berdistribusi normal atau memiliki hubungan non-linear antara variabel respons dan prediktor. GLM sering digunakan untuk data kategorik, count, atau data proporsi, dan mencakup banyak metode populer seperti regresi logistik dan regresi Poisson.

GLM terdiri dari 3 komponen utama :

Distribusi (Random Component):
Menentukan distribusi probabilitas dari variabel respons (Y). Tidak harus normal.
Fungsi Link (Systematic Component):

Fungsi yang menghubungkan nilai harapan E(Y) dengan kombinasi linear dari prediktor.
Prediktor Linear (Linear Predictor):
Kombinasi linear dari variabel bebas:

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

8.1 Exponential Family

Distribusi disebut termasuk dalam exponential family jika bisa dituliskan dalam bentuk umum berikut:

\[ f(y; \theta) = h(y) \exp\left[\eta(\theta) \cdot T(y) - A(\theta)\right] \]

Di mana:

$y$ = data atau variabel respons
$\theta$ = parameter
$\eta(\theta)$ = natural parameter (parameter kanonik)
$T(y)$ = sufficient statistic
$A(\theta)$ = fungsi normalisasi (cumulant function)
$h(y)$ = fungsi dasar yang hanya tergantung pada $y$

Contoh distribusi yang termasuk eksponential family :

Normal (Gaussian)
Bernoulli
Binomial
Poisson
Gamma (dengan parameter tertentu)
Exponential

Contoh pembuktian distribusi Poisson :

Fungsi massa probabilitas (pmf) dari Poisson:

\[ P(Y = y) = \frac{\mu^y e^{-\mu}}{y!}, \quad y = 0, 1, 2, \dots \]

Kita ubah bentuknya :

\[ P(Y = y) = \frac{1}{y!} \exp \left[ y \log \mu - \mu \right] \]

Bandingkan dengan bentuk umum:

8.2 Model Regresi Logistik

Regresi logistik adalah metode statistik yang digunakan untuk memodelkan hubungan antara satu atau lebih variabel prediktor (X) dengan variabel respon biner (Y) — yaitu variabel yang hanya memiliki dua kemungkinan hasil, seperti:

Ya/Tidak
Sakit/Sehat
Lulus/Tidak Lulus

Model ini memprediksi probabilitas suatu kejadian menggunakan fungsi logit:

\[ \log\left(\frac{p}{1 - p}\right) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k \]

Dengan:

p = probabilitas kejadian (misal: “ya”)
Sisi kiri = log odds

Klasifikasi Biner (Respon Biner)

Medis:
Memprediksi apakah pasien terjangkit penyakit (1) atau tidak (0) berdasarkan usia, tekanan darah, dan kolesterol.
Pendidikan:
Menentukan apakah siswa akan lulus (1) atau tidak (0) dari suatu pelatihan berdasarkan nilai awal dan kehadiran.
Marketing:
Menentukan apakah konsumen akan membeli produk (1) atau tidak (0) berdasarkan jenis kelamin, umur, dan penghasilan.
HR & Rekrutmen:
Prediksi apakah seorang kandidat akan diterima (1) atau tidak (0) berdasarkan nilai psikotes dan pengalaman kerja.

Keunggulan Utama Regresi Logistik:
1. Menghasilkan probabilitas (0–1), bukan hanya klasifikasi.
  → Cocok untuk pengambilan keputusan berbasis risiko.
2. Tidak mengharuskan hubungan linear antara X dan Y, seperti pada regresi linear biasa.
3. Cocok untuk variabel prediktor campuran: numerik, kategorik, atau gabungan keduanya.
4. Interpretasi parameter mudah melalui odds ratio.
5. Mudah diperluas ke regresi logistik multikategori atau regresi logistik ordinal.
Fungsi Sigmoid

Fungsi sigmoid adalah fungsi matematis yang mengubah nilai apa pun di garis bilangan real menjadi angka antara 0 dan 1. Ini membuatnya sangat cocok untuk model klasifikasi biner, karena hasil akhirnya bisa ditafsirkan sebagai probabilitas dari sebuah kejadian.

\[ \sigma(z) = \frac{1}{1 + e^{-z}} \]
- $z$ adalah nilai input (bisa berupa hasil prediksi linier: $\beta_0 + \beta_1 x_1 + \cdots$)
- $e$ adalah bilangan Euler (sekitar 2.718)
- Outputnya akan selalu berada dalam rentang 0 sampai 1

Estimasi Parameter

data = read.csv(file.choose(), header= TRUE, sep = ",")
data

##     balance housing
## 1      2143       1
## 2        29       1
## 3         2       1
## 4      1506       1
## 5         1       0
## 6       231       1
## 7       447       1
## 8         2       1
## 9       121       1
## 10      593       1
## 11      270       1
## 12      390       1
## 13        6       1
## 14       71       1
## 15      162       1
## 16      229       1
## 17       13       1
## 18       52       1
## 19       60       1
## 20        0       1
## 21      723       1
## 22      779       1
## 23       23       1
## 24       50       1
## 25        0       1
## 26     -372       1
## 27      255       1
## 28      113       1
## 29     -246       1
## 30      265       1
## 31      839       0
## 32      378       1
## 33       39       1
## 34        0       1
## 35    10635       1
## 36       63       1
## 37       -7       1
## 38       -3       0
## 39      506       1
## 40        0       1
## 41     2586       1
## 42       49       1
## 43      104       1
## 44      529       1
## 45       96       1
## 46     -171       1
## 47     -364       1
## 48        0       1
## 49        0       0
## 50        0       1
## 51     1291       1
## 52     -244       1
## 53        0       1
## 54      -76       1
## 55     -103       1
## 56      243       0
## 57      424       1
## 58      306       1
## 59       24       1
## 60      179       1
## 61        0       1
## 62      989       1
## 63      249       1
## 64      790       1
## 65      154       1
## 66     6530       1
## 67      100       0
## 68       59       1
## 69     1205       1
## 70    12223       1
## 71     5935       1
## 72       25       1
## 73      282       1
## 74       23       1
## 75     1937       1
## 76      384       1
## 77      582       0
## 78       91       0
## 79        0       1
## 80        1       1
## 81      206       1
## 82      164       0
## 83      690       1
## 84     2343       1
## 85      137       1
## 86      173       1
## 87       45       0
## 88     1270       1
## 89       16       1
## 90      486       1
## 91       50       0
## 92      152       1
## 93      290       1
## 94       54       1
## 95      -37       0
## 96      101       1
## 97      383       0
## 98       81       1
## 99        0       1
## 100     229       1
## 101    -674       1
## 102      90       0
## 103     128       1
## 104     179       1
## 105       0       1
## 106      54       1
## 107     151       1
## 108      61       0
## 109      30       1
## 110     523       1
## 111      31       1
## 112      79       0
## 113     -34       1
## 114     448       1
## 115      81       1
## 116     144       1
## 117     351       1
## 118     -67       1
## 119     262       0
## 120       0       0
## 121      56       1
## 122      26       1
## 123       3       1
## 124      41       1
## 125       7       0
## 126     105       1
## 127     818       1
## 128     -16       1
## 129       0       1
## 130    2476       1
## 131    1185       0
## 132     217       1
## 133    1685       1
## 134     802       1
## 135       0       1
## 136      94       1
## 137       0       1
## 138       0       0
## 139     517       1
## 140     265       1
## 141     947       1
## 142       3       1
## 143      42       0
## 144      37       1
## 145      57       1
## 146      22       1
## 147       8       1
## 148     293       1
## 149       3       1
## 150     348       1

y = data$housing
x = data$balance

modellog = glm(y ~ x, data = data, family = binomial)
modellog

## 
## Call:  glm(formula = y ~ x, family = binomial, data = data)
## 
## Coefficients:
## (Intercept)            x  
##   1.6125493    0.0005059  
## 
## Degrees of Freedom: 149 Total (i.e. Null);  148 Residual
## Null Deviance:       125.1 
## Residual Deviance: 123   AIC: 127

Visualisasi

library(ggplot2)

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

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

8.3 Model Regresi Poisson

Distribusi Poisson memiliki fungsi probabilitas (pmf):

\[ P(Y = y) = \frac{e^{-\lambda} \lambda^y}{y!}, \quad y = 0, 1, 2, \ldots \]

dapati ditulis dalam bentuk exponential family :

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

dengan bentuk umu exponential family :

\[ f(y; \theta) = \exp \left( y\theta - b(\theta) + c(y) \right) \]

maka :

$\theta$ (canonical parameter) = $\log(\lambda)$
$b(\theta) = \lambda = \exp(\theta)$
$c(y) = -\log(y!)$

Fungsi Link Canonical

Dalam GLM, kita menghubungkan mean dari Y (λ) dengan prediktor linier melalui fungsi link.

Untuk fungsi link canonical adalah :

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

\[ \log(\mu_i) = \beta_0 + \beta_1 x_{1i} + \dots + \beta_k x_{ki} \]atau :

\[ \mu_i = \exp(\beta_0 + \beta_1 x_{1i} + \dots + \beta_k x_{ki}) \]Estimasi Regresi Poisson

poisson_model = glm(y ~ x, data = data, family = poisson)
summary(poisson_model)

## 
## Call:
## glm(formula = y ~ x, family = poisson, data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -1.697e-01  9.313e-02  -1.822   0.0685 .
## x            2.090e-05  5.265e-05   0.397   0.6913  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 40.603  on 149  degrees of freedom
## Residual deviance: 40.454  on 148  degrees of freedom
## AIC: 300.45
## 
## Number of Fisher Scoring iterations: 4

Visualisasi

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

9. Inferensi GLM

Dalam model GLM, kita sering ingin mengetahui seberapa baik parameter yang kita estimasi menggambarkan kondisi sebenarnya. Untuk itu, kita perlu memahami dua hal penting: rata-rata dari estimasi (ekspektasi) dan sebaran estimasi (varians). Pengetahuan ini jadi dasar untuk uji statistik seperti uji Wald, Likelihood Ratio Test, atau dalam membuat interval kepercayaan.

1. Ekspektasi Estimator

Ekspektasi menjelaskan apakah suatu estimator memiliki bias atau tidak. Sebuah estimator dikatakan tak bias jika nilai harapannya sama dengan nilai parameter sebenarnya. Dalam notasi matematika:

\[ E(\hat{\beta}) = \beta \]2. Varians Estimator

Varians menggambarkan tingkat ketepatan dari estimasi parameter. Jika varians kecil, berarti estimasi tersebut konsisten dan tidak mudah berubah-ubah. Varians dari Beta dapat dihitung secara mendekati dengan rumus berikut:

Distribusi Asimptotik

Untuk sampel besar, penyebaran dari estimasi parameter Beta dapat didekati dengan distribusi normal sebagai berikut:

\[ E(\hat{\beta}) = \beta \]Distribusi ini menjadi dasar untuk berbagai inferensi seperti :

Menguji apakah suatu koefisien signifikan atau tidak (uji Wald),
Membuat interval kepercayaan untuk Beta
Menghitung p-value dalam uji hipotesis.

model = glm(y ~ x, family = poisson)
summary(model)

## 
## Call:
## glm(formula = y ~ x, family = poisson)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -1.697e-01  9.313e-02  -1.822   0.0685 .
## x            2.090e-05  5.265e-05   0.397   0.6913  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 40.603  on 149  degrees of freedom
## Residual deviance: 40.454  on 148  degrees of freedom
## AIC: 300.45
## 
## Number of Fisher Scoring iterations: 4

9.1 Mencari Ekspetasi dan Varians dalam GLM

Ekspetasi

Jika suatu distribusi termasuk keluarga eksponensial:

\[ \log f(y; \theta) = y \cdot \theta - b(\theta) + c(y) \]Maka ekspektasi dari Y diperoleh dari turunan pertama fungsi log-partisi:

\[ E(Y) = b'(\theta) \]Varians

Varians dari Y diperoleh dari turunan kedua fungsi log-partisi:

\[ \mathrm{Var}(Y) = b''(\theta) \]Jika ada parameter dispersi φ, bentuknya jadi:

9.2 Metode Penaksiran Parameter

Dalam Generalized Linear Model (GLM), salah satu tujuan utama adalah melakukan penaksiran terhadap parameter model, yaitu koefisien regresi yang menghubungkan variabel prediktor dengan rata-rata respon. Penaksiran ini dilakukan berdasarkan pendekatan Maximum Likelihood Estimation (MLE).

Metode MLE digunakan untuk mencari nilai parameter yang memaksimalkan fungsi likelihood, yaitu peluang munculnya data yang diamati berdasarkan distribusi yang diasumsikan.

Langkah-langkah Penaksiran MLE dalam GLM

Tentukan bentuk fungsi likelihood dari data, berdasarkan distribusi keluarga eksponensial.
Ambil log-likelihood untuk menyederhanakan perhitungan:

\[ \ell(\beta) = \log L(\beta) \]
Turunkan log-likelihood terhadap parameter β:

\[ \frac{\partial \ell(\beta)}{\partial \beta} = 0 \]

Ini disebut sebagai score equation.
Gunakan metode iteratif, karena score equation umumnya tidak bisa diselesaikan secara analitik. Salah satu metode paling populer adalah:
- Iteratively Reweighted Least Squares (IRLS)

9.3 Diagnostik Model GLM

Setelah membangun model GLM, langkah penting selanjutnya adalah melakukan diagnostik model untuk mengevaluasi sejauh mana model tersebut cocok dengan data. Tujuannya adalah untuk mendeteksi potensi masalah seperti pencilan (outlier), pengaruh besar (influential points), atau kesalahan spesifikasi model.

1. Deviance

Deviance adalah ukuran seberapa baik model memodelkan data dibandingkan dengan model yang sepenuhnya cocok (model saturasi).

\[ D = 2 \left[ \ell(\hat{y}_{\text{saturated}}) - \ell(\hat{y}_{\text{model}}) \right] \]

$\ell(\cdot)$ adalah fungsi log-likelihood
Semakin kecil nilai deviance, semakin baik model

2. Pearson Residuals

Pearson residuals adalah selisih antara nilai observasi dan nilai prediksi yang dinormalisasi.

\[ r_i^{(P)} = \frac{y_i - \hat{\mu}_i}{\sqrt{V(\hat{\mu}_i)}} \]

$y_i$ = nilai aktual
$\hat{\mu}_i$ = nilai harapan (fitted)
$V(\hat{\mu}_i)$ = varians dari $y_i$

3. Deviance Residual

Deviance residual untuk observasi ke- i adalah:

\[ r_i^{(D)} = \text{sign}(y_i - \hat{\mu}_i) \cdot \sqrt{d_i} \]

9.4 Detail Metode Estimasi dan Inferensi Regresi Logistik

Untuk memperoleh nilai parameter β pada regresi logistik, kita perlu memaksimalkan fungsi log-likelihood dari model. Salah satu pendekatan numerik yang umum digunakan adalah metode Newton-Raphson, karena model tidak memiliki solusi eksplisit.

Fungsi Log-Likelihood pada Regresi Logistik

Model regresi logistik memodelkan probabilitas keberhasilan sebagai berikut:

\[ \pi_i = \frac{1}{1 + \exp(-\mathbf{x}_i^\top \boldsymbol{\beta})} \]Untuk n data observasi, fungsi log-likelihood ditulis sebagai :

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

Langkah-Langkah Estimasi Newton-Raphson pada Regresi Logistik

Gradient (Turunan Pertama) / Score Function

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

Turunan Kedua (Hessian Matrix)

\[ H(\boldsymbol{\beta}) = -\mathbf{X}^\top \mathbf{W} \mathbf{X} \]

dengan:

\[ \mathbf{W} = \text{diag}(\pi_i (1 - \pi_i)) \]

Formula Iterasi Newton-Raphson

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

Keterangan:

$\boldsymbol{\pi}^{(t)}$: vektor probabilitas pada iterasi ke-$t$
Iterasi dilakukan hingga nilai parameter $\boldsymbol{\beta}$ konvergen.

# Variabel prediktor dan respons
x <- as.numeric(data$balance)
y <- as.numeric(data$housing)
n <- length(y)

# Tambah intercept
X <- cbind(1, x)  # Matriks desain dengan intercept
beta <- matrix(c(0, 0), ncol = 1)  # Inisialisasi koefisien

# Fungsi log-likelihood
loglik <- function(beta, X, y) {
  eta <- X %*% beta
  mu <- exp(eta)
  sum(y * eta - mu - log(factorial(y)))
}

# Newton-Raphson Iterasi
tol <- 1e-6
max_iter <- 100
for (i in 1:max_iter) {
  eta <- X %*% beta
  mu <- exp(eta)
  
  # Gradient (turunan pertama)
  grad <- t(X) %*% (y - mu)
  
  # Hessian (turunan kedua)
  W <- diag(as.vector(mu))
  H <- -t(X) %*% W %*% X
  
  # Update parameter
  beta_new <- beta - solve(H) %*% grad
  
  # Cek konvergensi
  if (sum(abs(beta_new - beta)) < tol) {
    cat("Konvergen pada iterasi ke-", i, "\n")
    break
  }
  
  beta <- beta_new
}

## Konvergen pada iterasi ke- 4

# Output akhir
cat("Koefisien hasil Newton-Raphson:\n")

## Koefisien hasil Newton-Raphson:

print(beta)

##            [,1]
##   -1.696736e-01
## x  2.090356e-05

Interpretasi :

Intercept (-0.16967): Saat balance = 0, maka rata-rata kejadian (housing) diperkirakan sekitar exp(-0.16967) ≈ 0.844.

Koefisien balance (0.0000209): Setiap peningkatan 1 unit pada balance meningkatkan ekspektasi kejadian sebesar exp(0.0000209) ≈ 1.000021, atau naik sekitar 0.0021%.

Uji Wald

Mengukur signifikansi statistik dari suatu parameter dengan membandingkan nilai estimasi parameter terhadap simpangan bakunya.

\[ W = \left( \frac{\hat{\beta} - \beta_0}{\text{SE}(\hat{\beta})} \right)^2, \quad W \sim \chi^2_1 \]

# Matriks Hessian terakhir dari Newton-Raphson
eta <- X %*% beta
mu <- exp(eta)
W <- diag(as.vector(mu))
H <- -t(X) %*% W %*% X

# Kovarians dari beta (inverse dari -Hessian)
cov_beta <- solve(-H)

# Standard error (akar diagonal kovarians)
se_beta <- sqrt(diag(cov_beta))

# Wald statistic
wald_stat <- (beta[,1] / se_beta)^2

# p-value dari distribusi chi-square dengan df=1
p_values <- 1 - pchisq(wald_stat, df = 1)

# Buat tabel hasil uji Wald
wald_result <- data.frame(
  Estimate = beta[,1],
  Std_Error = se_beta,
  Wald_Stat = wald_stat,
  p_value = p_values
)
rownames(wald_result) <- c("Intercept", "balance")

# Tampilkan hasil
wald_result

##                Estimate    Std_Error Wald_Stat    p_value
## Intercept -1.696736e-01 9.313020e-02 3.3193065 0.06847071
## balance    2.090356e-05 5.264906e-05 0.1576375 0.69134101

Interpretasi :

Berdasarkan uji Wald, tidak ada bukti kuat bahwa baik intercept maupun variabel balance berpengaruh secara signifikan terhadap variabel dependen pada tingkat signifikansi 5%.

Uji Likelihood Ratio

Menguji apakah model dengan parameter tambahan secara signifikan lebih baik dari model yang lebih sederhana.

\[ LR = -2 \left[ \ell(\text{model terbatas}) - \ell(\text{model penuh}) \right], \quad LR \sim \chi^2_k \]

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

# Model penuh (full model): dengan prediktor balance
model_full <- glm(housing ~ balance, data = data, family = binomial)

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

## Analysis of Deviance Table
## 
## Model 1: housing ~ 1
## Model 2: housing ~ balance
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       149     125.06                     
## 2       148     123.01  1   2.0552   0.1517

Interpretasi :

Variabel balance tidak berpengaruh signifikan terhadap variabel housing berdasarkan uji Likelihood Ratio (p-value = 0.1517).

9.5 Detail Metode Estimasi dan Inferensi Regresi Poisson

Distribusi Poisson

Distribusi Poisson memiliki fungsi probabilitas:

\[ P(Y = y) = \frac{\lambda^y e^{-\lambda}}{y!}, \quad y = 0, 1, 2, \dots \]

Di mana $\lambda$ adalah rata-rata dan juga varians dari distribusi Poisson.

Model Regresi Poisson

Dalam regresi Poisson, kita memodelkan rata-rata $\lambda_i$ sebagai:

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

Sehingga:

\[ Y_i \sim \text{Poisson}(\lambda_i) \]

Fungsi Log-Likelihood Poisson

Fungsi log-likelihood untuk distribusi Poisson didefinisikan sebagai:

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

Dengan:

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

Sehingga, fungsi log-likelihood regresi Poisson dapat ditulis sebagai:

\[ \ell(\beta) = \sum_{i=1}^n \left[ y_i (\mathbf{x}_i^\top \beta) - \exp(\mathbf{x}_i^\top \beta) - \log(y_i!) \right] \]

Estimasi dilakukan dengan metode iterasi (IRLS).

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

Tahap 1: Definisikan Model Regresi Poisson

Model regresi Poisson didefinisikan sebagai:

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

Tahap 2: Log-Likelihood yang Dimaksimalkan

Fungsi log-likelihood untuk distribusi Poisson adalah:

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

Tahap 3: Formulasi Iteratif (IRLS)

Langkah iteratif untuk memperbarui parameter $\boldsymbol{\beta}$ adalah:

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

Dengan:

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

dan

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

Pengujian Hipotesis

Uji Wald

# Koefisien dan standar error
coef_val <- coef(model)[2]
se_val <- summary(model)$coefficients[2, 2]

wald_z <- coef_val / se_val
wald_chisq <- wald_z^2
p_value <- 1 - pchisq(wald_chisq, df = 1)

cat("Z:", wald_z, "\nChi-Square:", wald_chisq, "\np-value:", p_value)

## Z: 0.397036 
## Chi-Square: 0.1576376 
## p-value: 0.6913409

Interpretasi :

Karena p-value > 0.05, maka terima H0 sehingga tidak ada pengaruh signifikan variabel balance terhadap kemungkinan memiliki housing, pada tingkat signifikansi 5%

Uji Likelihood Ratio (Chi-Square)

# Koefisien dan standar error
coef_val <- coef(model)[2]
se_val <- summary(model)$coefficients[2, 2]

wald_z <- coef_val / se_val
wald_chisq <- wald_z^2
p_value <- 1 - pchisq(wald_chisq, df = 1)

cat("Z:", wald_z, "\nChi-Square:", wald_chisq, "\np-value:", p_value)

## Z: 0.397036 
## Chi-Square: 0.1576376 
## p-value: 0.6913409

Interpretasi :

Karena p-value > 0.05, maka terima H0 sehingga tidak ada pengaruh signifikan variabel balance dalam menjelaskan variabel housing dalam model, pada taraf signifikansi 5%

Evaluasi Model (AIC & BIC)

AIC(model)

## [1] 300.4545

BIC(model)

## [1] 306.4757

10. Regresi Logistik dengan Prediktor Normal, Ordinal, dan Rasio

Regresi logistik adalah metode statistik untuk memodelkan hubungan antara satu atau lebih variabel prediktor dengan variabel dependen kategori biner (misalnya: sukses/gagal, ya/tidak). Model ini memperkirakan probabilitas terjadinya suatu kejadian berdasarkan nilai prediktor.

Skala prediktor yang umum digunakan dalam regresi logistik:

Nominal Variabel kategorik tanpa urutan atau tingkatan. Contohnya: jenis kelamin (laki-laki/perempuan), warna favorit (merah, biru, hijau). Biasanya dikonversi ke variabel dummy agar bisa digunakan dalam model.
Ordinal Variabel kategorik yang memiliki urutan atau tingkatan, tetapi jarak antar kategori tidak harus sama. Contohnya: tingkat pendidikan (SMA, Sarjana, Magister), skala kepuasan (rendah, sedang, tinggi).
- Rasio Variabel numerik kontinu yang memiliki nol mutlak dan jarak antar nilai bermakna. Contohnya: usia, pendapatan, berat badan. Variabel ini dapat langsung digunakan dalam model tanpa perlu konversi.

10.1 Simulasi Data

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

## # A tibble: 6 × 4
##   resign gender job_level workload
##    <int> <chr>  <chr>        <dbl>
## 1      0 Male   Junior        50.1
## 2      1 Male   Junior        49.6
## 3      1 Male   Senior        45  
## 4      1 Male   Manager       45.7
## 5      0 Female Junior        41.4
## 6      1 Female Manager       44

10.2 Eksplorasi Data

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:vcdExtra':
## 
##     summarise

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

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

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

sim_data %>%
  group_by(resign) %>%
  summarise(
    Jumlah = n(),
    Rata2_Workload = mean(workload)
  )

## # A tibble: 2 × 3
##   resign Jumlah Rata2_Workload
##    <int>  <int>          <dbl>
## 1      0     65           44.4
## 2      1    435           45.0

10.3 Perlakuan Variabel Ordinal

10.3.1 Treat sebagai Nominal (Dummy)

 sim_data_nominal <- sim_data %>%
 mutate(
 job_level = factor(job_level, levels = c("Junior", "Senior", "Manager"))
 )
 model_nominal <- glm(resign ~ gender + job_level + workload, data = sim_data_nominal, family = binomial)
 summary(model_nominal)

## 
## Call:
## glm(formula = resign ~ gender + job_level + workload, family = binomial, 
##     data = sim_data_nominal)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       0.91009    1.17380   0.775  0.43814   
## genderMale       -0.39802    0.27166  -1.465  0.14287   
## job_levelSenior   0.69283    0.32715   2.118  0.03419 * 
## job_levelManager  1.36235    0.48849   2.789  0.00529 **
## workload          0.01865    0.02613   0.714  0.47541   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 386.39  on 499  degrees of freedom
## Residual deviance: 371.64  on 495  degrees of freedom
## AIC: 381.64
## 
## Number of Fisher Scoring iterations: 5

Interpretasi Koefisien:

Intercept (0.91009):
Saat semua variabel prediktor bernilai nol (misalnya gender = perempuan, job level = junior, dan workload = 0), maka log-odds dari kejadian yang dimodelkan adalah 0.91. Namun, nilai ini tidak signifikan secara statistik (p = 0.43814).
genderMale (-0.39802):
Menjadi laki-laki menurunkan log-odds dari kejadian dibanding perempuan sebesar 0.398. Tetapi efek ini tidak signifikan (p = 0.14287), jadi tidak ada bukti kuat bahwa gender berpengaruh.
job_levelSenior (0.69283):
Menjadi pegawai senior (dibandingkan dengan level referensi, yaitu junior) meningkatkan log-odds dari kejadian sebesar 0.69. Efek ini signifikan secara statistik (p = 0.03419), menunjukkan bahwa level senior berpengaruh terhadap probabilitas kejadian.
job_levelManager (1.36235):
Menjadi manager meningkatkan log-odds dari kejadian sebesar 1.36 dibandingkan pegawai junior. Efek ini sangat signifikan (p = 0.00529), menunjukkan bahwa menjadi manager secara jelas meningkatkan peluang terjadinya kejadian.
workload (0.01865):
Setiap peningkatan satu unit pada workload sedikit meningkatkan log-odds dari kejadian sebesar 0.0186, namun efek ini tidak signifikan (p = 0.47541), sehingga tidak ada bukti kuat bahwa workload memengaruhi hasil.

Evaluasi Model:

Null deviance = 386.39 dan residual deviance = 371.64 menunjukkan sedikit penurunan, artinya penambahan prediktor memberikan sedikit perbaikan model.
AIC = 381.64, digunakan untuk membandingkan model lain — semakin rendah, semakin baik.
Iterasi Fisher Scoring = 5, menandakan bahwa model konvergen dengan baik dan tidak ada masalah dalam proses estimasi.

Kesimpulan:

Hanya variabel job level (senior dan manager) yang secara signifikan memengaruhi probabilitas terjadinya kejadian.
Variabel gender dan workload tidak menunjukkan pengaruh yang signifikan.
Model masih dapat ditingkatkan, terutama karena penurunan deviance dari model null ke model akhir tidak terlalu besar.

10.3.2 Treat sebagai Rasio

 sim_data_numeric <- sim_data %>%
 mutate(
 job_level_numeric = case_when(
 job_level == "Junior" ~ 1,
 job_level == "Senior" ~ 2,
 job_level == "Manager" ~ 3
 )
 )
 model_numeric <- glm(resign ~ gender + job_level_numeric + workload, data = sim_data_numeric, family = binomial)
 summary(model_numeric)

## 
## Call:
## glm(formula = resign ~ gender + job_level_numeric + workload, 
##     family = binomial, data = sim_data_numeric)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        0.22478    1.19660   0.188  0.85100   
## genderMale        -0.39769    0.27144  -1.465  0.14289   
## job_level_numeric  0.68504    0.21121   3.243  0.00118 **
## workload           0.01868    0.02611   0.715  0.47447   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 386.39  on 499  degrees of freedom
## Residual deviance: 371.64  on 496  degrees of freedom
## AIC: 379.64
## 
## Number of Fisher Scoring iterations: 5

Interpretasi Koefisien:

(Intercept) = 0.22478
Ketika semua variabel prediktor bernilai nol (misalnya: gender = perempuan, job_level_numeric = 0, workload = 0), log-odds seseorang untuk resign adalah 0.22478. Namun, ini tidak signifikan (p = 0.85100), jadi tidak memberikan informasi kuat.
genderMale = -0.39769
Menjadi laki-laki menurunkan log-odds untuk resign sebesar 0.397 dibanding perempuan. Namun, efek ini tidak signifikan (p = 0.14289), artinya tidak ada bukti kuat bahwa gender memengaruhi keputusan resign.
job_level_numeric = 0.68504
Setiap kenaikan satu tingkat jabatan (misalnya dari junior ke senior, atau senior ke manager) meningkatkan log-odds seseorang untuk resign sebesar 0.685. Efek ini signifikan secara statistik (p = 0.00118), menunjukkan bahwa semakin tinggi level jabatan, semakin besar kecenderungan untuk resign.
workload = 0.01868
Setiap peningkatan satu unit workload hanya sedikit meningkatkan log-odds resign sebesar 0.0187. Efek ini tidak signifikan (p = 0.47447), jadi tidak ada bukti kuat bahwa workload berpengaruh terhadap resign.

Evaluasi Model:

Null deviance = 386.39
Ini adalah deviance model tanpa prediktor (hanya intercept).
Residual deviance = 371.64
Ini adalah deviance model dengan semua prediktor. Penurunan ini menunjukkan adanya perbaikan, meskipun tidak drastis.
AIC = 379.64
Digunakan untuk membandingkan model lain. Lebih rendah = lebih baik.
Jumlah iterasi = 5
Model berhasil konvergen dalam 5 iterasi, artinya estimasi stabil.

Kesimpulan :

Satu-satunya variabel yang signifikan adalah job_level_numeric, yang menunjukkan bahwa makin tinggi jabatan seseorang, makin besar peluangnya untuk resign.
Gender dan workload tidak signifikan, sehingga belum ada cukup bukti bahwa dua variabel ini berpengaruh terhadap keputusan resign.
Model ini bisa digunakan sebagai dasar analisis, tapi masih mungkin diperbaiki dengan eksplorasi variabel lain atau interaksi antar variabel.

Visualisasi Prediksi

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

Interpretasi:

Pola Linear Positif untuk Semua Job Level
Untuk setiap job_level, seiring meningkatnya workload, probabilitas resign (asumsinya, berdasarkan model resign ~ ...) sedikit meningkat. Artinya, beban kerja yang lebih tinggi sedikit meningkatkan kemungkinan resign, namun efeknya tidak terlalu curam.
Perbedaan Probabilitas antar Job Level:
- Manager (biru) memiliki probabilitas tertinggi untuk resign, bahkan pada beban kerja yang rendah.
- Senior (hijau) ada di tengah.
- Junior (merah) memiliki probabilitas paling rendah untuk resign, konsisten di seluruh rentang workload.
- Ini menunjukkan bahwa model memperkirakan semakin tinggi level jabatan, semakin besar kecenderungan untuk resign. Hal ini konsisten dengan hasil regresi sebelumnya, di mana variabel job_level_numeric signifikan positif.
Workload Tidak Terlalu Berpengaruh:
Kenaikan probabilitas karena workload tampak hampir datar (kemiringan lemah), sesuai dengan hasil regresi sebelumnya di mana workload tidak signifikan secara statistik. Jadi model tetap memperkirakan sedikit peningkatan, tetapi tidak tajam.
Job Level Sebagai Faktor Dominan:
Model secara jelas membedakan probabilitas berdasarkan level jabatan. Terlihat ada tiga “lapisan” titik terpisah (Junior → Senior → Manager), menunjukkan bahwa job_level merupakan prediktor yang jauh lebih kuat dibanding workload.

Kesimpulan:

Faktor job level jauh lebih penting dalam memprediksi probabilitas resign dibanding workload.
Workload hanya memberi efek kecil, meskipun naik sedikit saat workload meningkat.
Model menangkap pola ordinal dalam job level (Junior < Senior < Manager), meskipun variabel job_level dalam model ini diperlakukan sebagai nominal.

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

Interpretasi Utama:

Pola Peningkatan Probabilitas Seiring Workload

-   Untuk setiap level jabatan, peningkatan `workload` (beban kerja) disertai dengan **kenaikan probabilitas untuk resign**.

-   Ini menunjukkan bahwa dalam model ini, **workload dipandang sebagai variabel numerik yang berkontribusi positif**, meskipun dari hasil regresi sebelumnya, efeknya **tidak signifikan secara statistik**.

Level Jabatan Berpengaruh Konsisten & Berurutan:

-    **Probabilitas resign meningkat secara berjenjang dari Junior → Senior → Manager**.

-    Ini sangat sesuai dengan penggunaan **job_level sebagai numerik (ordinal)** karena:

    -    Probabilitas diprediksi **secara teratur meningkat** dari level 1 ke 3.

    -   **Jarak antar grup lebih simetris dan teratur** dibanding model sebelumnya (yang memperlakukan job_level sebagai nominal).

-   Ini memperkuat temuan dari regresi bahwa **job_level_numeric signifikan**, dan model mengenali hubungan ordinalnya dengan baik.

Manager Tetap Paling Rentan Resign:

-   Sama seperti pada grafik sebelumnya, **Manager (level 3)** memiliki probabilitas resign tertinggi, bahkan pada beban kerja rendah.

-    **Junior (level 1)** tetap paling kecil kemungkinannya untuk resign.

Kesimpulan:

Memperlakukan job level sebagai numerik (ordinal) menghasilkan pola prediksi yang lebih halus dan konsisten.
Workload masih berkontribusi menaikkan probabilitas resign, meskipun efeknya tidak dominan.
Model ini lebih mencerminkan hubungan hirarkis antara jabatan dan kemungkinan resign, menjadikannya lebih cocok digunakan dibanding model dengan job_level sebagai nominal.

Ringkasan Koefisien Model

library(knitr)
library(kableExtra)

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

## 
## Attaching package: 'kableExtra'

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

library(broom)

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

# Ringkasan model nominal
tidy(model_nominal) %>%
  kable(format = "latex", 
        booktabs = TRUE, 
        caption = "Ringkasan Koefisien Model dengan Ordinal sebagai Nominal") %>%
  kable_styling(latex_options = c("hold_position", "striped"))

11. Pemilihan Model Regresi Logistik dan Evaluasi

11.1 Membangun Model Regresi Logistik : Pendekatan Confirmatory dan Exploratory

1. Pendekatan Confirmatory (Konfirmatori)

Tujuan:

Menguji hipotesis atau teori yang telah ada sebelumnya. Model dibangun berdasarkan landasan teori atau temuan riset sebelumnya.

Ciri-ciri:

Variabel bebas dipilih berdasarkan teori atau kajian ilmiah sebelumnya.
Peneliti sudah memiliki dugaan hubungan antara variabel bebas dan variabel dependen.
Lebih fokus pada pengujian signifikansi hubungan tertentu.
Tidak terlalu banyak eksplorasi data.
Sering digunakan dalam penelitian akademik atau pengujian teori.

Contoh Penggunaan:

Seorang peneliti ingin menguji apakah dukungan sosial dan status pekerjaan mempengaruhi kemungkinan seseorang mengalami depresi, berdasarkan teori psikologi yang sudah ada. Maka, model logistik dibangun hanya dengan variabel yang sesuai teori tersebut.

2. Pendekatan Explanatory (Eksplanatori)

Tujuan:

Menjelaskan atau menemukan variabel-variabel yang paling berpengaruh terhadap variabel dependen. Model dibangun dengan pendekatan yang lebih eksploratif.

Ciri-ciri:

Variabel bebas dapat dipilih berdasarkan eksplorasi data (misalnya melalui uji korelasi, seleksi variabel otomatis).
Bertujuan untuk menemukan pola atau prediktor yang signifikan.
Sering menggunakan metode seleksi variabel (forward, backward, stepwise).
Digunakan dalam data mining atau penelitian awal (eksploratif).
Fokus pada akurasi model, bukan hanya pengujian teori.

Contoh Penggunaan:

Sebuah perusahaan e-commerce ingin memprediksi peluang pelanggan melakukan pembelian ulang. Mereka memasukkan banyak variabel seperti frekuensi login, jumlah belanja terakhir, jenis kelamin, umur, lokasi, dsb., lalu menggunakan stepwise logistic regression untuk menemukan prediktor paling kuat, tanpa hipotesis awal.

Perbandingan Singkat:

Aspek	Confirmatory	Explanatory
Tujuan utama	Menguji teori/hipotesis	Menemukan pola/prediktor baru
Dasar pemilihan variabel	Teori atau penelitian sebelumnya	Eksplorasi data/statistik
Sifat analisis	Deduktif	Induktif
Digunakan dalam	Penelitian akademik	Praktik bisnis, data mining
Teknik seleksi variabel	Manual (berdasarkan teori)	Otomatis (stepwise, lasso, dsb.)

Simulasi Data

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

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

## Loading required package: lattice

## 
## Attaching package: 'lattice'

## The following object is masked from 'package:gnm':
## 
##     barley

## 
## Attaching package: 'caret'

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

library(pROC)

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

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

## 
## Attaching package: 'pROC'

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

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

##   y          x1 x2         x3
## 1 0  0.07716537  0 -0.6810714
## 2 0  0.37431557  1 -0.3265440
## 3 1 -0.33481390  1 -0.9318383
## 4 0 -0.94988897  0 -0.3925825
## 5 0 -0.94021008  0  0.1528529
## 6 1  1.18962312  1 -0.6612526

Pemilihan Model

Model Full

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

## 
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.8237     0.2462  -3.346  0.00082 ***
## x1            1.0367     0.2107   4.919 8.69e-07 ***
## x2           -0.3256     0.3432  -0.949  0.34268    
## x3            0.5357     0.1952   2.745  0.00606 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 244.35  on 199  degrees of freedom
## Residual deviance: 205.81  on 196  degrees of freedom
## AIC: 213.81
## 
## Number of Fisher Scoring iterations: 4

Penjelasan masing-masing:

Intercept (-0.8237): Ketika semua prediktor bernilai nol, maka log-odds (log peluang) dari y = 1 adalah -0.8237. Dalam konteks probabilitas, ini berarti probabilitas dasar cukup rendah.
x1 (1.0367): Koefisien signifikan secara statistik (p < 0.001). Ini berarti setiap kenaikan 1 unit pada x1 akan meningkatkan log-odds dari y = 1 sebesar 1.0367.
Dalam bentuk odds ratio:

\[ e^{1.0367} \approx 2.82 \]

Artinya, setiap kenaikan satu unit x1 meningkatkan peluang kejadian (odds) y = 1 sekitar 2.8 kali.
x2 (-0.3256): Tidak signifikan secara statistik (p = 0.343). Artinya, tidak cukup bukti bahwa x2 memengaruhi y.
x3 (0.5357): Signifikan (p = 0.00606). Setiap kenaikan satu unit pada x3 meningkatkan log-odds y = 1 sebesar 0.5357.
Odds ratio-nya:

\[ e^{0.5357} \approx 1.71 \]

Jadi, peluang y = 1 meningkat sekitar 71% untuk setiap kenaikan satu unit x3.

Goodness of Fit:

Null deviance = 244.35 (model tanpa prediktor)
Residual deviance = 205.81 (model dengan x1, x2, x3)
Penurunan deviance sebesar 38.54 menunjukkan bahwa model dengan prediktor jauh lebih baik dibandingkan tanpa prediktor.
AIC = 213.81: Digunakan untuk membandingkan model. Semakin kecil, semakin baik.

Kesimpulan Umum:

Model ini signifikan secara keseluruhan (terlihat dari deviance yang turun).
x1 dan x3 adalah prediktor yang signifikan dan berkontribusi pada model.
x2 tidak signifikan, dan bisa dipertimbangkan untuk dieliminasi dalam model yang lebih parsimonious.
Model dapat digunakan untuk memprediksi probabilitas kejadian y = 1.

11.2 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)
 132

## [1] 132

AIC(model_full, step_forward, step_backward, step_both)

##               df      AIC
## model_full     4 213.8055
## step_forward   3 212.7091
## step_backward  3 212.7091
## step_both      3 212.7091

11.3 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")

11.4 Pseudo R-Squared

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

##   CoxSnell Nagelkerke   McFadden 
##  0.1715370  0.2432183  0.1540302

11.5 Tabel Klasifikasi dan Evaluasi

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

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 127  39
##          1  13  21
##                                           
##                Accuracy : 0.74            
##                  95% CI : (0.6734, 0.7993)
##     No Information Rate : 0.7             
##     P-Value [Acc > NIR] : 0.1227748       
##                                           
##                   Kappa : 0.2935          
##                                           
##  Mcnemar's Test P-Value : 0.0005265       
##                                           
##             Sensitivity : 0.3500          
##             Specificity : 0.9071          
##          Pos Pred Value : 0.6176          
##          Neg Pred Value : 0.7651          
##              Prevalence : 0.3000          
##          Detection Rate : 0.1050          
##    Detection Prevalence : 0.1700          
##       Balanced Accuracy : 0.6286          
##                                           
##        'Positive' Class : 1               
##

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

## Sensitivity Specificity 
##   0.3500000   0.9071429

11.6 Metode Perbandingan Model dalam Regresi Logistik

library(MASS)
 library(broom)
 library(DescTools)
 set.seed(69)

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

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

##     Model      AIC Deviance
## 1 Model 1 312.0390 308.0390
## 2 Model 2 299.3297 293.3297
## 3 Model 3 274.3247 266.3247

11.7 Likelihood-Ratio Test

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

## Analysis of Deviance Table
## 
## Model 1: y ~ x1
## Model 2: y ~ x1 + x2
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       298     308.04                          
## 2       297     293.33  1   14.709 0.0001254 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

## Analysis of Deviance Table
## 
## Model 1: y ~ x1 + x2
## Model 2: y ~ x1 + x2 + x3
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       297     293.33                          
## 2       296     266.32  1   27.005 2.029e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

11.8 Prinsip Parsimony

Rumus:

Rumus AIC :

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

Log(L) : log-likelihood dari model (semakin besar, model makin cocok).

k: jumlah parameter dalam model (termasuk intercept).

AIC menyeimbangkan kecocokan model dan kompleksitas model.

Rumus Deviance :

\[ \text{Deviance} = -2 \cdot \left[ \log(L_{\text{model}}) - \log(L_{\text{saturated}}) \right] \]

log(Lmodel): log-likelihood dari model yang kita uji.

log(Lsaturated): log-likelihood dari model saturasi (model dengan kecocokan sempurna).

Rumus Likelihood-Ratio

\[ \text{LR Statistic} = -2 \cdot \left[ \log(L_{\text{restricted}}) - \log(L_{\text{full}}) \right] \]

Lrestricted: likelihood dari model lebih kecil (misalnya hanya intercept).

Lfull: likelihood dari model lengkap (dengan semua prediktor).

Uji statistik ini mengikuti distribusi chi-kuadrat dengan derajat bebas = jumlah parameter tambahan di model penuh.

11.9 Evaluasi Tabel Klasifikasi dan Akurasi Model

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

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 195  43
##          1  19  43
##                                          
##                Accuracy : 0.7933         
##                  95% CI : (0.743, 0.8377)
##     No Information Rate : 0.7133         
##     P-Value [Acc > NIR] : 0.001011       
##                                          
##                   Kappa : 0.4487         
##                                          
##  Mcnemar's Test P-Value : 0.003489       
##                                          
##             Sensitivity : 0.5000         
##             Specificity : 0.9112         
##          Pos Pred Value : 0.6935         
##          Neg Pred Value : 0.8193         
##              Prevalence : 0.2867         
##          Detection Rate : 0.1433         
##    Detection Prevalence : 0.2067         
##       Balanced Accuracy : 0.7056         
##                                          
##        'Positive' Class : 1              
##

11.9.1 Sensitivitas dan Spesifisitas

Sensitifitas

Definisi:

Ukuran kemampuan model dalam mendeteksi kasus positif dengan benar.

\[ \text{Sensitivitas} = \frac{\text{True Positive (TP)}}{\text{True Positive (TP)} + \text{False Negative (FN)}} \]

Spesitifitas

Definisi:

Ukuran kemampuan model dalam mendeteksi kasus negatif dengan benar.

\[ \text{Spesifisitas} = \frac{\text{True Negative (TN)}}{\text{True Negative (TN)} + \text{False Positive (FP)}} \]

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

## Sensitivity Specificity 
##    0.500000    0.911215

11.10 Detail ROC Penjelasan Kurva ROC (Receiver Operating Characteristic)

1. Definisi Kurva ROC

ROC curve adalah grafik yang menggambarkan kinerja model klasifikasi biner dengan memplot:

\[ \text{True Positive Rate (TPR)} = \text{Sensitivitas} \]

\[ \text{False Positive Rate (TPR)} = \text{1 - Spesifitas} \]

Sumbu X = False Positive Rate (FPR)
Sumbu Y = True Positive Rate (TPR)

2. Cut-off dan Pergerakan Kurva

Cut-off (Threshold):

Dalam regresi logistik, probabilitas default dipotong di 0.5 untuk menentukan positif/negatif.
Namun, ROC menggambarkan performa model pada semua nilai cut-off, dari 0 hingga 1.

Pergerakan Kurva:

Saat cut-off diturunkan (misalnya dari 0.9 ke 0.1), model lebih sering memprediksi positif, sehingga:
- TPR (sensitivitas) meningkat
- FPR juga meningkat
Setiap titik pada kurva ROC menunjukkan kinerja model pada suatu cut-off tertentu.

3. Kurva ROC Ideal

Kurva ideal: melewati sudut kiri atas (TPR = 1, FPR = 0)
Ini berarti:

-    Semua kasus positif terdeteksi (TPR = 1)

-    Tidak ada kasus negatif yang salah klasifikasi (FPR = 0)

Semakin dekat ke pojok kiri atas, semakin baik modelnya.

4. Interpretasi Luas Area di Bawah Kurva (AUC - Area Under Curve)

Nilai AUC:

AUC	Interpretasi Model
0.90 – 1.00	Sangat baik (excellent)
0.80 – 0.90	Baik (good)
0.70 – 0.80	Cukup (fair)
0.60 – 0.70	Buruk (poor)
0.50 – 0.60	Gagal, hampir acak
= 0.50	Tidak lebih baik dari tebak-tebakan (random)
< 0.50	Model lebih buruk dari acak (bermasalah)

Makna AUC:

Probabilitas bahwa model akan memberikan nilai skor yang lebih tinggi untuk kasus positif dibanding kasus negatif secara acak.

5. Kegunaan Kurva ROC

Evaluasi performa model secara umum

Tidak bergantung pada cut-off tertentu.
Cocok ketika distribusi kelas tidak seimbang.

Membandingkan model

Model dengan AUC lebih tinggi lebih baik dalam membedakan antara kelas.

Memilih threshold terbaik

Bisa digunakan bersama dengan kurva lain (misal: kurva Youden Index atau Precision-Recall) untuk menentukan cut-off optimal.

6.Visualisasi dalam R

 library(pROC)
 set.seed(69)
 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)

7.Simulasi Pemilihan Threshold Optimal

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

##    Threshold Sensitivity Specificity
## 1       0.10  0.94339623   0.4693878
## 2       0.15  0.92452830   0.5918367
## 3       0.20  0.86792453   0.6938776
## 4       0.25  0.79245283   0.7619048
## 5       0.30  0.75471698   0.8163265
## 6       0.35  0.67924528   0.8503401
## 7       0.40  0.62264151   0.8571429
## 8       0.45  0.60377358   0.9047619
## 9       0.50  0.43396226   0.9251701
## 10      0.55  0.43396226   0.9251701
## 11      0.60  0.37735849   0.9387755
## 12      0.65  0.30188679   0.9591837
## 13      0.70  0.22641509   0.9727891
## 14      0.75  0.16981132   0.9863946
## 15      0.80  0.16981132   0.9931973
## 16      0.85  0.09433962   0.9931973
## 17      0.90  0.09433962   1.0000000

8.Catatan

ROC cocok saat proporsi kelas seimbang

Untuk data yang tidak seimbang, dapat menggunakan precision-recall curve

11.11 Precision-Recall Curve (PR Curve)

1. Definisi PR Curve

PR Curve (Precision-Recall Curve) adalah grafik yang digunakan untuk mengevaluasi kinerja model klasifikasi biner, terutama saat data imbalance (tidak seimbang).

Sumbu-sumbu PR Curve:

X-axis: Recall (Sensitivitas)

$\frac{TP}{TP + FP}$
Y-axis: Precision

$\frac{TP}{TP + FP}$

2. Interpretasi PR Curve

PR Curve menunjukkan trade-off antara precision dan recall untuk berbagai nilai threshold (cut-off).
Saat threshold berubah, precision dan recall akan berubah, dan titik-titik tersebut membentuk kurva.X

Makna bentuk kurva:

Kurva mendekati pojok kanan atas (precision & recall tinggi) → model sangat bagus.
Kurva mendekati garis dasar (precision rendah walau recall tinggi) → banyak false positives, model buruk.

3. Area Under PR Curve (AUPRC / AP)

AUPRC (Area Under the Precision-Recall Curve) atau Average Precision (AP) adalah nilai ringkasan kinerja model.
Semakin besar AUPRC, semakin baik model dalam mendeteksi kelas positif secara tepat.

Interpretasi AUPRC:

Nilai antara 0 hingga 1
Nilai tinggi (> 0.8) menunjukkan precision tetap tinggi walau recall meningkat
Lebih sensitif daripada AUC ROC saat data tidak seimbang

4. PR Curve vs ROC Curve

Aspek	PR Curve	ROC Curve
Sumbu X	Recall (TPR)	False Positive Rate (FPR)
Sumbu Y	Precision	True Positive Rate (TPR)
Fokus utama	Kinerja di kelas positif	Kinerja di semua kelas
Cocok untuk	Data imbalanced (kelas minoritas)	Data seimbang
Kurva bagus	Dekat pojok kanan atas	Dekat pojok kiri atas
Nilai area	AUPRC (lebih rendah dari AUC saat imbalance)	AUC (sering terlihat bagus meski model lemah di minoritas)
Contoh kasus	Deteksi penipuan, kanker, spam	Umum, data seimbang

5.Visualisasi PR Curve di R

library(PRROC)

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

## Loading required package: rlang

 set.seed(69)
 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)

11.12 Pseudo R-squared pada Regresi Logistik

Tujuan

Tujuan Pseudo R-squared pada regresi logistik adalah untuk memberikan ukuran goodness-of-fit (kebaikan model) — yaitu, seberapa baik model menjelaskan variabilitas data — meskipun tidak sebanding langsung dengan R^2 pada regresi linear.

Simulasi Data

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

Model Logistik dan Null Model

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

Likelihood dan Rumus

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

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

dengan :

$L_0$ : likelihood model null (tanpa prediktor)
$L_M$ : likelihood model null (tanpa prediktor)

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

##   R2_Cox_Snell R2_McFadden
## 1    0.2669425   0.2591474

Perhitungan Otomatis dengan Package Tambahan

Menggunakan pscl

library(pscl)

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

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

pR2(model)

## fitting null model for pseudo-r2

##          llh      llhNull           G2     McFadden         r2ML         r2CU 
## -133.1623342 -179.7420077   93.1593470    0.2591474    0.2669425    0.3822818

Menggunakann rcompanion

library(rcompanion)

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

 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.259147
## Cox and Snell (ML)                   0.266943
## Nagelkerke (Cragg and Uhler)         0.382282
## 
## $Likelihood.ratio.test
##  Df.diff LogLik.diff  Chisq    p.value
##       -3      -46.58 93.159 4.5904e-20
## 
## $Number.of.observations
##           
## Model: 300
## Null:  300
## 
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
## 
## $Warnings
## [1] "None"

Menggunakan DescTools

library(DescTools)
 PseudoR2(model, which = "all")

##        McFadden     McFaddenAdj        CoxSnell      Nagelkerke   AldrichNelson 
##       0.2591474       0.2368933       0.2669425       0.3822818       0.2369506 
## VeallZimmermann           Efron McKelveyZavoina            Tjur             AIC 
##       0.4346929       0.2939840       0.4377139       0.2951128     274.3246684 
##             BIC          logLik         logLik0              G2 
##     289.1397983    -133.1623342    -179.7420077      93.1593470

Kesimpulan Sederhana:

Model memiliki performa cukup baik: McFadden > 0.2 dan Nagelkerke > 0.3.
Performa model lebih baik dibanding model null secara signifikan (lihat G²)
AIC & BIC bisa digunakan untuk membandingkan dengan model lain.

12. Apa itu Distribusi Multinomial

Distribusi multinomial adalah perluasan dari distribusi binomial ke kasus dengan lebih dari dua kategori hasil.

Digunakan ketika suatu percobaan (eksperimen) menghasilkan salah satu dari k kategori yang saling eksklusif.
Dilakukan sebanyak n kali (percobaan independen).
Probabilitas setiap kategori tetap konstan dalam tiap percobaan.

$P(X_1 = x_1, \dots, X_k = x_k) = \frac{n!}{x_1! \cdots x_k!} \cdot p_1^{x_1} \cdots p_k^{x_k}$

dengan :

$n$ : jumlah total percobaan

$k$ : jumlah kategori (kelas hasil)

$x_i$ : jumlah kejadian pada kategori ke-$i$, dengan $\sum_{i=1}^k x_i = n$

$p_i$ : probabilitas kategori ke-$i$, dengan $\sum_{i=1}^k p_i = 1$
$X_i$ : variabel acak jumlah kemunculan kategori ke-$i$

12.1 Studi Kasus

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

Studi Kasus:

Sebuah toko makanan ringan mencatat preferensi 100 pelanggan terhadap 3 jenis snack:

Keripik (A): probabilitas disukai $p_1 = 0.4$
Biskuit (B): probabilitas disukai $p_2 = 0.35$
Permen (C): probabilitas disukai $p_3 = 0.25$
Dari 100 pelanggan, ternyata:
45 memilih Keripik ($x_1 = 45$)
30 memilih Biskuit ($x_2 = 30$)
25 memilih Permen ($x_3 = 25$)

Pertanyaan:

Berapa probabilitas (peluang) bahwa distribusi pilihan pelanggan seperti di atas dapat terjadi?

Rumus Distribusi Multinomial

\[ P(x_1 = 45, x_2 = 30, x_3 = 25) = \frac{100!}{45! \cdot 30! \cdot 25!} \cdot 0.4^{45} \cdot 0.35^{30} \cdot 0.25^{25} \]

Perhitungan Manual di R

# Data input
n <- 100
x1 <- 45
x2 <- 30
x3 <- 25
p1 <- 0.4
p2 <- 0.35
p3 <- 0.25

# Hitung komponen log faktorial agar tidak overflow
log_n_fact <- lfactorial(n)
log_x1_fact <- lfactorial(x1)
log_x2_fact <- lfactorial(x2)
log_x3_fact <- lfactorial(x3)

# Hitung log probabilitas
log_prob <- log_n_fact - (log_x1_fact + log_x2_fact + log_x3_fact) +
  x1 * log(p1) + x2 * log(p2) + x3 * log(p3)

# Ambil eksponensial untuk dapatkan probabilitas sebenarnya
probabilitas <- exp(log_prob)
print(probabilitas)

## [1] 0.004376503

Jadi, peluang terjadinya 45 orang memilih Keripik, 30 Biskuit, dan 25 Permen adalah sekitar 4.84%.

12.2.1 Baseline-category logit model

Baseline-category logit model (juga dikenal sebagai multinomial logit model) digunakan untuk menganalisis variabel respon kategorik dengan lebih dari dua kategori tanpa urutan (nominal).

Model ini membandingkan setiap kategori terhadap satu kategori referensi (baseline) menggunakan logit (log odds).

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

dengan :

$\pi_j$

Probabilitas bahwa suatu pengamatan termasuk ke dalam kategori $j$

$\pi_c$

Probabilitas baseline (kategori referensi)

$\alpha_j$

Intersep logit untuk kategori $j$

$\beta_j$

Koefisien regresi logit kategori $j$

$x$

Variabel prediktor (bisa skalar atau vektor)

$j = 1, \dots, c-1$

Semua kategori kecuali kategori dasar ($c$)

Estimasi Parameter

Estimasi parameter pada Baseline-Category Logit Model (juga dikenal sebagai Multinomial Logit Model) dilakukan dengan pendekatan Log-Likelihood

\[ \log L(\boldsymbol{\theta}) = \sum_{i=1}^{n} \sum_{j=1}^{c} y_{ij} \log(\pi_{ij}) \]

12.3 Contoh Kasus

Sebuah universitas ingin mengetahui faktor-faktor yang memengaruhi pemilihan moda transportasi mahasiswa saat pergi ke kampus.
Moda transportasi yang diteliti: Motor, Mobil, atau Transportasi Umum.

Universitas melakukan survei terhadap 200 mahasiswa dan mengumpulkan data berikut:

Transport: Moda transportasi utama yang digunakan (Motor, Mobil, Umum)
Income: Pendapatan bulanan mahasiswa (dalam juta rupiah)
Gender: Jenis kelamin (Male, Female)
Distance: Jarak dari rumah ke kampus (dalam km)

Tujuannya: Mengetahui bagaimana pendapatan, jenis kelamin, dan jarak tempat tinggal memengaruhi pilihan moda transportasi mahasiswa.

12.4 Simulasi Data

set.seed(69)

n <- 200

# Variabel prediktor
Gender <- sample(c("Male", "Female"), n, replace = TRUE)
Income <- round(rnorm(n, mean = 3, sd = 1), 1)         # dalam juta
Distance <- round(rnorm(n, mean = 7, sd = 3), 1)        # dalam km

# Simulasi pilihan transportasi
Transport <- mapply(function(g, inc, dist) {
  if (inc >= 4) {
    sample(c("Mobil", "Motor", "Umum"), 1, prob = c(0.6, 0.3, 0.1))
  } else if (dist > 10) {
    sample(c("Umum", "Motor", "Mobil"), 1, prob = c(0.5, 0.3, 0.2))
  } else {
    sample(c("Motor", "Umum", "Mobil"), 1, prob = c(0.5, 0.4, 0.1))
  }
}, Gender, Income, Distance)

# Buat data frame
df <- data.frame(
  Transport = factor(Transport),
  Income = Income,
  Gender = factor(Gender),
  Distance = Distance
)

# Jadikan "Motor" sebagai baseline
df$Transport <- relevel(df$Transport, ref = "Motor")

# Tampilkan data awal
head(df)

##   Transport Income Gender Distance
## 1     Mobil    5.2 Female      5.0
## 2     Mobil    4.6   Male      6.0
## 3     Mobil    3.9   Male      4.2
## 4     Motor    3.8 Female      5.8
## 5      Umum    1.9 Female      7.5
## 6      Umum    3.1 Female      5.0

12.5 Estimasi Model

library(nnet)

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

# Estimasi model
model <- multinom(Transport ~ Income + Gender + Distance, data = df)

## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 187.873104
## final  value 187.440548 
## converged

# Ringkasan model
summary(model)

## Call:
## multinom(formula = Transport ~ Income + Gender + Distance, data = df)
## 
## Coefficients:
##       (Intercept)     Income   GenderMale   Distance
## Mobil  -5.9089243  1.1891218 -0.006775961 0.12627122
## Umum   -0.3059968 -0.1284675  0.146577936 0.06473189
## 
## Std. Errors:
##       (Intercept)    Income GenderMale   Distance
## Mobil   1.1376922 0.2544691  0.4384373 0.07678836
## Umum    0.6739959 0.1745558  0.3150413 0.05554210
## 
## Residual Deviance: 374.8811 
## AIC: 390.8811

12.6 Nilai p-value dan Interpretasi

# Z-value dan p-value
z <- summary(model)$coefficients / summary(model)$standard.errors
pval <- 2 * (1 - pnorm(abs(z)))
round(pval, 4)

##       (Intercept) Income GenderMale Distance
## Mobil      0.0000 0.0000     0.9877   0.1001
## Umum       0.6498 0.4617     0.6417   0.2438

Interpretasi :

Kategori “Mobil”

Income: 0.0000 → Sangat signifikan
GenderMale: 0.9877 → Tidak signifikan
Distance: 0.1001 → Tidak signifikan (mendekati signifikan, tapi masih > 0.05)

Kategori “Umum”

Income: 0.4617 → Tidak signifikan
GenderMale: 0.6417 → Tidak signifikan
Distance: 0.2438 → Tidak signifikan

Untuk kategori Mobil, hanya Income yang signifikan dalam mempengaruhi pilihan ke kategori ini.

Untuk kategori Umum, tidak ada variabel yang signifikan pada level 5%.

12.7 Prediksi dan Validasi

df$Predicted <- predict(model)
 table(Predicted = df$Predicted, Actual = df$Transport)

##          Actual
## Predicted Motor Mobil Umum
##     Motor    69    13   56
##     Mobil     6    17    5
##     Umum     15     5   14

12.8 Kesimpulan

Memprediksi variabel dependen kategori dengan lebih dari dua kelas
→ Digunakan ketika variabel respon (Y) bersifat kategorik non-ordinal dengan >2 kategori, misalnya: jenis transportasi (Mobil, Motor, Umum), jenis pekerjaan, atau preferensi produk.

Menganalisis pengaruh variabel independen terhadap peluang memilih suatu kategori
→ Menunjukkan bagaimana variabel numerik atau kategorik (seperti pendapatan, jenis kelamin, jarak) memengaruhi probabilitas seseorang memilih salah satu kategori dari variabel respon.

Membandingkan kategori terhadap kategori referensi (baseline)
→ Model secara otomatis membandingkan setiap kategori terhadap satu kategori acuan, sehingga memungkinkan interpretasi relatif antar pilihan. Misalnya: seberapa besar pengaruh pendapatan dalam memilih “Mobil” dibanding “Motor”

12.8.1 Contoh Kasus 2 di R

Menggunakan dataset Iris

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

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

summary(model)

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

Interpretasi Koefisien

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

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

Prediksi dan Visualisasi

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

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

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

13 Regresi Logistik Ordinal

Regresi logistik ordinal digunakan ketika variabel dependen (respon) memiliki kategori yang berurutan secara logis, tapi jaraknya tidak bisa diukur secara pasti.

13.1 Konsep Cumulative Logit Model

\[ \log\left( \frac{P(Y \leq j)}{P(Y > j)} \right) = \theta_j - \boldsymbol{\beta}^\top x, \quad \text{untuk } j = 1, 2, \ldots, J - 1 \]

dengan :

Variabel respon ordinal: Y ∈ {1, 2, …, J}

$x$ : vektor variabel prediktor

$β$: koefisien regresi

$θ$ : ambang (threshold) untuk kategori ke-jjj

13.2 Interpreasi Koefisien

Setiap koefisien $\beta_k$ mewakili pengaruh variabel prediktor $x_k$ terhadap log-odds kumulatif (yaitu log peluang untuk berada di kategori lebih rendah atau sama dibanding lebih tinggi).

Jika $\beta_k$ > 0
→ Peluang berada di kategori yang lebih rendah meningkat,
→ atau sebaliknya, peluang berada di kategori lebih tinggi menurun saat $x_k$ naik.

Jika $\beta_k$ < 0
→ Peluang berada di kategori yang lebih tinggi meningkat,
→ atau sebaliknya, peluang berada di kategori lebih rendah menurun saat $x_k$ naik.

13.3 Contoh Data Kepuasan Penggunaan Layanan Transportasi

# Simulasi Data
set.seed(123)
n <- 200

# Kecepatan layanan (semakin tinggi semakin cepat)
speed <- round(runif(n, 1, 10))

# Kepuasan pelanggan (semakin tinggi speed → cenderung lebih puas)
satisfaction <- cut(
  5 + 0.5 * speed + rnorm(n),
  breaks = c(-Inf, 5.5, 7.5, Inf),
  labels = c("Tidak Puas", "Cukup", "Puas"),
  ordered_result = TRUE
)

# Gabung jadi data frame
df <- data.frame(satisfaction, speed)
head(df)

##   satisfaction speed
## 1        Cukup     4
## 2         Puas     8
## 3        Cukup     5
## 4         Puas     9
## 5         Puas     9
## 6   Tidak Puas     1

13.4 Estimasi Model Ordinal

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

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

13.5 Nilai P-Value

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

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

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

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

13.6 Prediksi Probabilitas

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

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

13.7 Goodness-of-Fit dan Proportional Odds

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

13.8 Alternatif Model Ordinal

Selain cumulative logit, model ordinal lainnya:

Adjacent-category logit
Continuation-ratio (sequential) logit

13.9 Kesimpulan

Digunakan untuk variabel respon yang bersifat kategorik dan berurutan (ordinal)

Contohnya: tingkat kepuasan (Tidak Puas, Cukup, Puas), tingkat pendidikan, atau level risiko.
Menganalisis bagaimana variabel prediktor memengaruhi peluang respon berada pada kategori yang lebih tinggi atau lebih rendah

Model memperkirakan log odds kumulatif dari kategori respon dan mengasumsikan efek prediktor konsisten di seluruh level (asumsi proportional odds).
Koefisien regresi diinterpretasikan sebagai perubahan log odds (atau odds ratio) terhadap peluang berada di kategori yang lebih tinggi
- Koefisien positif → meningkatkan peluang ke kategori lebih tinggi
- Koefisien negatif → menurunkan peluang ke kategori lebih tinggi

13.10 Asumsi Paralelisme dalam Regresi Logistik Ordinal

Asumsi ini menyatakan bahwa:

Efek dari setiap variabel prediktor adalah konstan di seluruh batas (cutpoint) kategori respon.

Dalam model:

\[ \log\left( \frac{P(Y \leq j)}{P(Y > j)} \right) = \theta_j - \boldsymbol{\beta}^\top x \]

Koefisien $\beta$ sama untuk semua $j$ (kategori).

Konsekuensi Pelanggaran

Jika asumsi ini tidak terpenuhi, maka:

Interpretasi koefisien menjadi tidak valid
Model bisa menjadi bias atau menyesatkan
Pilihan model ordinal tidak lagi tepat → disarankan pakai generalized ordered logit atau multinomial logit (kalau urutan tidak terlalu penting)

Pengujian Asumsi Paralelisme

Beberapa metode umum:

a. Brant Test (di R: brant() dari package brant)

Menguji apakah koefisien signifikan berbeda antar cutpoint
Hasil:

-   **p-value \< 0.05** → asumsi dilanggar (tidak paralel)

-    **p-value \> 0.05** → asumsi terpenuhi

b. Likelihood Ratio Test (LRT)

Bandingkan model ordinal (dengan asumsi) vs model non-paralel (tanpa asumsi)
Signifikan → asumsi tidak valid

Kesimpulan

Asumsi paralelisme menyatakan bahwa efek setiap prediktor terhadap peluang kumulatif antar kategori adalah konstan.
Artinya, koefisien regresi untuk semua batas kategori (cutpoint) diasumsikan sama.

Jika asumsi ini dilanggar, maka model ordinal menjadi tidak valid dan interpretasi koefisien bisa menyesatkan.
Dalam kasus seperti ini, disarankan menggunakan model alternatif seperti generalized ordered logit atau multinomial logit.

Asumsi ini dapat diuji dengan uji Brant atau likelihood ratio test (LRT).
Jika hasil pengujian menunjukkan p-value < 0.05, maka asumsi paralelisme tidak terpenuhi.

14 Log Linear Model

Log-linear model adalah model statistik yang digunakan untuk menganalisis hubungan antar variabel kategorik yang disajikan dalam bentuk tabel kontingensi. Model ini tidak membedakan antara variabel dependen dan independen, karena semua variabel dianggap setara. Tujuan utama dari log-linear model adalah untuk memahami struktur asosiasi atau interaksi antar variabel dalam data kategorik, serta mengevaluasi apakah variabel-variabel tersebut saling bebas (independen) atau tidak.

Dalam log-linear model, kita memodelkan logaritma dari frekuensi yang diharapkan dalam setiap sel tabel kontingensi sebagai fungsi dari efek utama (main effects) dan efek interaksi. Misalnya, untuk dua variabel kategorik $A$ dan $B$, model independensi ditulis sebagai:

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

di mana $\mu_{ij}$ adalah nilai harapan (expected count) untuk sel ke- $i$, $j$ .Jika interaksi dimasukkan, maka ditambah komponen $\lambda^{AB}_{ij}$, sehingga menjadi model saturasi.

Log-linear model dapat bersifat saturated, artinya model memasukkan semua efek interaksi (2 arah, 3 arah, dst.) dan akan memodelkan data dengan sempurna (deviance = 0). Namun, dalam praktiknya, kita lebih sering menggunakan model yang direduksi (reduced model) untuk menghindari overfitting dan untuk menguji apakah pengaruh tertentu (seperti interaksi) signifikan. Perbandingan antara model saturasi dan model reduksi dilakukan melalui uji deviance atau likelihood ratio test.

Setiap parameter dalam log-linear model berhubungan dengan log-count dari data kategorik. Parameter efek utama $\lambda^A_i$ menunjukkan kontribusi level ke- $i$ dari variabel A terhadap frekuensi sel. Parameter interaksi, seperti $\lambda^{AB}_{ij}$, menunjukkan apakah terdapat ketergantungan antara kombinasi level $i$ dari A dan $j$ dari B. Jika interaksi signifikan, maka variabel A dan B tidak independen.

Log-linear model paling umum digunakan untuk menganalisis tabel kontingensi 2-dimensi, 3-dimensi, atau lebih. Misalnya, untuk data tiga variabel seperti Jenis Kelamin (G), Status Merokok (S), dan Penyakit Paru (D), kita bisa membentuk model log-linear seperti:

\[ \log(\mu_{ijk}) = \lambda + \lambda_i^G + \lambda_j^S + \lambda_k^D + \lambda_{ij}^{GS} + \lambda_{ik}^{GD} + \lambda_{jk}^{SD} \]

Untuk menguji apakah ada interaksi antara ketiga variabel, kita tambahkan $\lambda^{GSD}_{ijk}$ ke dalam model.

Berbeda dengan regresi logistik, di mana ada satu variabel sebagai outcome (biner atau multinomial), log-linear model tidak memposisikan variabel mana pun sebagai respon. Oleh karena itu, log-linear lebih cocok digunakan jika kita tertarik pada simetri hubungan antar variabel kategorik, bukan untuk prediksi. Namun, dalam kasus 2x2 atau 2xk tabel, log-linear model dan logistik biner dapat menghasilkan informasi yang setara.

Tujuan utama dari log-linear model adalah untuk mengevaluasi struktur ketergantungan antar variabel kategorik, menentukan model terbaik yang menjelaskan data tanpa overfitting, dan membantu dalam reduksi dimensi interaksi yang tidak signifikan. Model ini juga digunakan sebagai dasar dalam analisis lanjutan seperti model logit (untuk data multinomial) dan model hierarki dalam analisis tabel kontingensi kompleks.

Contoh :

data_glm <- data.frame(
  Lulus = c(1, 0, 1, 0),
  Pelatihan = factor(c("Ya", "Ya", "Tidak", "Tidak")),
  Frek = c(80, 20, 50, 50)
)

# Model logistik biner dengan frekuensi sebagai weights
model_logit <- glm(Lulus ~ Pelatihan, weights = Frek, family = binomial, data = data_glm)

# Tampilkan hasil
summary(model_logit)

## 
## Call:
## glm(formula = Lulus ~ Pelatihan, family = binomial, data = data_glm, 
##     weights = Frek)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.081e-16  2.000e-01    0.00        1    
## PelatihanYa  1.386e+00  3.202e-01    4.33 1.49e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 258.98  on 3  degrees of freedom
## Residual deviance: 238.71  on 2  degrees of freedom
## AIC: 242.71
## 
## Number of Fisher Scoring iterations: 5

Aspek	Tabel Kontingensi	Model Log-linear	Regresi Logistik
1. Tujuan Analisis	Mendeskripsikan frekuensi observasi untuk kombinasi kategori	Menjelaskan asosiasi dan interaksi antar variabel kategorik	Memprediksi probabilitas kejadian berdasarkan satu variabel respons dan kovariat
2. Struktur Model	Tidak menggunakan model matematis	Model log-frekuensi: semua variabel sebagai prediktor (tanpa variabel respon)	Model logit: satu variabel sebagai respons, sisanya sebagai prediktor
3. Jenis Variabel	Semua variabel kategorik	Semua variabel kategorik	Respons: kategorik (biner/multinomial); prediktor bisa kategorik atau numerik
4. Interpretasi Parameter	Tidak ada parameter	Parameter menyatakan efek utama dan interaksi log-frekuensi	Parameter menyatakan log-odds atau odds ratio dari variabel terhadap kejadian
5. Aplikasi Umum	Analisis deskriptif sederhana	Uji independensi/interaksi antar kategori dalam data cross-classified	Prediksi/pemodelan outcome biner atau multinomial berdasarkan faktor

14.1 Tabel Kontingensi dan Model Loglinier

matrix(c(80, 20, 50, 50), nrow=2,
 dimnames = list(Pelatihan = c("Ya","Tidak"),
 Lulus = c("Ya", "Tidak")))

##          Lulus
## Pelatihan Ya Tidak
##     Ya    80    50
##     Tidak 20    50

Model Log Linear untuk tabel I x J dapat ditulis :

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

dengan :

$λ$: intercept (konstanta)

$\lambda^A$: efek utama dari level ke- $i$ pada variabel $A$

$\lambda^{B}_{ij}$: efek utama dari level ke- $j$ pada variabel $B$

$\lambda^{AB}_{ij}$: efek interaksi antara level ke- $i$ dari $A$ dan ke- $j$ dari $B$

14.2 Model Saturated

Model yang memasukkan semua efek utama dan seluruh interaksi antar variabel dalam data, sehingga model mampu memprediksi data dengan sempurna.

Contoh :

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

##         Serangan Ya Tidak
## Obat                     
## Timolol          35    65
## Placebo          45    55

 model_saturated <- loglm(~ Obat * Serangan, data = data)
 summary(model_saturated)

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

14.3 Model Independent

Model ini mengasumsikan tidak adanya interaksi

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

 model_indep <- loglm(~ Obat + Serangan, data = data)
 summary(model_indep)

## Formula:
## ~Obat + Serangan
## attr(,"variables")
## list(Obat, Serangan)
## attr(,"factors")
##          Obat Serangan
## Obat        1        0
## Serangan    0        1
## attr(,"term.labels")
## [1] "Obat"     "Serangan"
## attr(,"order")
## [1] 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 0
## attr(,".Environment")
## <environment: R_GlobalEnv>
## 
## Statistics:
##                       X^2 df  P(> X^2)
## Likelihood Ratio 2.087576  1 0.1485015
## Pearson          2.083333  1 0.1489147

14.4 Odds Ratio dan Interpretasi

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

14.5 Estimasi Parameter

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

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

## [1] -0.4183685

14.6 Model Lebih Sederhana dan Perbandingan Model

 anova(model_indep, model_saturated)

## LR tests for hierarchical log-linear models
## 
## Model 1:
##  ~Obat + Serangan 
## Model 2:
##  ~Obat * Serangan 
## 
##           Deviance df Delta(Dev) Delta(df) P(> Delta(Dev)
## Model 1   2.087576  1                                    
## Model 2   0.000000  0   2.087576         1         0.1485
## Saturated 0.000000  0   0.000000         0         1.0000

14.7 Studi Kasus: Hubungan antara Kepuasan Kerja dan Keinginan Pindah Kerja

# Buat matriks data
data_survey <- matrix(c(45, 75,
                        80, 520,
                        25, 578),
                      nrow = 3, byrow = TRUE,
                      dimnames = list(Kepuasan = c("Tidak Puas", "Cukup Puas", "Sangat Puas"),
                                      PindahKerja = c("Ya", "Tidak")))

# Tampilkan tabel
ftable(data_survey)

##             PindahKerja  Ya Tidak
## Kepuasan                         
## Tidak Puas               45    75
## Cukup Puas               80   520
## Sangat Puas              25   578

15 MODEL LOG LINEAR 2 ARAH

15.1 Pengertian Model Log Linear pada Tabel Kontigensi

Model log-linear adalah model statistik yang digunakan untuk menganalisis hubungan (asosiasi/interaksi) antara dua atau lebih variabel kategorik dalam tabel kontingensi. Model ini bekerja dengan mengambil logaritma dari frekuensi harapan pada setiap sel dalam tabel, lalu mengekspresikannya sebagai kombinasi linier dari parameter yang merepresentasikan efek utama dan interaksi antar variabel.

15.2 Perbedaan Utama antara Model Log-Linear dan Model Regresi Logistik

Aspek	Log-Linear Model	Regresi Logistik
Tujuan	Menjelaskan asosiasi/interaksi antar variabel kategorik	Memprediksi probabilitas kejadian suatu kategori respon
Struktur Data	Tabel kontingensi (frekuensi sel)	Data individual (baris = unit responden/objek)
Variabel Respons	Tidak ada, semua variabel simetris	Ada satu variabel respon (biner atau multinomial)
Interpretasi Parameter	Efek terhadap log frekuensi sel	Efek terhadap log odds (rasio peluang)
Contoh Aplikasi	Uji independensi atau interaksi antar kategori	Prediksi probabilitas suatu hasil berdasarkan faktor

15.3 Analisis Data Tabel Kontingensi 2x2

# Membuat tabel kontingensi
tabel <- matrix(c(30, 10, 20, 40), 
                nrow = 2, byrow = FALSE,
                dimnames = list(Merokok = c("Ya", "Tidak"),
                                Kesehatan = c("Sakit", "Sehat")))

# Menampilkan tabel
tabel

##        Kesehatan
## Merokok Sakit Sehat
##   Ya       30    20
##   Tidak    10    40

15.4 Bentuk Model Log-Linear

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

dengan constraint sum-to-zero :

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

15.5 Estimasi Parameter Model (Manual, Sum-to-zero)

Langkah 1: Ambil nilai log frekuensi

Cell	Keterangan	Nilai	Log(Frek)
$\mu_{11}$	Merokok = Ya, Sakit	30	log(30)≈3.401
$\mu_{12}$	Merokok = Ya, Sehat	20	log(20)≈2.996
$\mu_{21}$	Merokok = Tidak, Sakit	10	log(10)≈2.303
$\mu_{22}$	Merokok = Tidak, Sehat	40	log(40)≈3.689

Langkah 2: Bentuk sistem persamaan

$\log(\mu_{22})$ == $λ$ (karena semua parameter lainnya nol di baseline)

$\lambda = \log(40) =3.689$

$\log(\mu_{12})$ = $\lambda$ + $\lambda_1^A$

$ _1^A$ = 2.996 - 3.689 = -0.693

$\log(\mu_{21})$ = $\lambda$ + $\lambda_1^B$

$\lambda_1^B$ = 2.303 - 3.689 = -1.386

$\log(\mu_{11})$ = $\lambda$ + $\lambda_1^A$ + $\lambda_1^B$ + $\lambda_{11}^{AB}$

$\lambda_{11}^{AB}$ = 3.401=3.689−0.693−1.386 =1.791

Hasil Estimasi Parameter:

$\lambda$ = 3.689
$\lambda_1^A$ = -0.693
$\lambda_1^B$ = -1.386
$\lambda_{11}^{AB}$ = 1.791

15.6 Hitung Odds Ratio dan Interval Kepercayaan

Langkah-langkah:

Log OR:

$\log(6)$ =  1.792

Standard Error (SE):

$$
SE = \sqrt{\frac{1}{30} + \frac{1}{20} + \frac{1}{10} + \frac{1}{40}} 
= \sqrt{0.0333 + 0.05 + 0.1 + 0.025} 
= \sqrt{0.2083} \approx 0.4564
$$

Interval di skala log:

$$
1.792 \pm 1.96 \times 0.4564 \Rightarrow 1.792 \pm 0.8946
$$

Batas bawah & atas di skala log:

-   Lower: 1.792 − 0.8946 = 0.8974

-   Upper: 1.792 + 0.8946 = 2.6866

Kembali ke skala OR (eksponensial):

\[ \text{CI } 95\% = \left( e^{0.8974},\, e^{2.6866} \right) \approx (2.45,\, 14.68) \]

15.7 Fitting Model Log-Linear dengan R

#Data2x2
 tabel<-matrix(c(30,20, 10, 40), nrow= 2,byrow=TRUE)
 colnames(tabel)<-c("Sakit", "Sehat")
 rownames(tabel)<-c("Ya","Tidak")
 tabel

##       Sakit Sehat
## Ya       30    20
## Tidak    10    40

data <-as.data.frame(as.table(tabel))
 colnames(data) <-c("Merokok", "Status", "Freq")
 data

##   Merokok Status Freq
## 1      Ya  Sakit   30
## 2   Tidak  Sakit   10
## 3      Ya  Sehat   20
## 4   Tidak  Sehat   40

# Model tanpa interaksi
 fit_no_inter <- glm(Freq ~ Merokok + Status, family = poisson, data = data)
 summary(fit_no_inter)

## 
## Call:
## glm(formula = Freq ~ Merokok + Status, family = poisson, data = data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.996e+00  1.871e-01  16.013   <2e-16 ***
## MerokokTidak 3.892e-10  2.000e-01   0.000    1.000    
## StatusSehat  4.055e-01  2.041e-01   1.986    0.047 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 21.288  on 3  degrees of freedom
## Residual deviance: 17.261  on 1  degrees of freedom
## AIC: 43.036
## 
## Number of Fisher Scoring iterations: 4

 # Model dengan interaksi
 fit_inter <- glm(Freq ~ Merokok * Status, family = poisson, data = data)
 summary(fit_inter)

## 
## Call:
## glm(formula = Freq ~ Merokok * Status, family = poisson, data = data)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.4012     0.1826  18.629  < 2e-16 ***
## MerokokTidak              -1.0986     0.3651  -3.009  0.00262 ** 
## StatusSehat               -0.4055     0.2887  -1.405  0.16015    
## MerokokTidak:StatusSehat   1.7918     0.4564   3.926 8.65e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.1288e+01  on 3  degrees of freedom
## Residual deviance: 3.9968e-15  on 0  degrees of freedom
## AIC: 27.775
## 
## Number of Fisher Scoring iterations: 3

15.8 Intepretasi Parameter

Parameter intercept merepresentasikan rata-rata logaritma dari frekuensi tiap sel.

Pengaruh dari variabel “Merokok” dan “Status” mencerminkan perbedaan nilai log frekuensi di antara kategori-kategori masing-masing variabel.

Adanya interaksi yang signifikan menunjukkan bahwa terdapat hubungan atau ketergantungan antara variabel Merokok dan Status kesehatan.

15.9 Analisis Data Tabel Kontingensi 2x3

Sebuah universitas ingin mengetahui hubungan antara jenis kelamin mahasiswa dan preferensi metode pembelajaran yang mereka sukai.

	Tatap Muka	Daring	Campuran
Laki-laki	15	18	7
Perempuan	22	20	10

15.10 Bentuk Model Log-Linear untuk Tabel 2x3

\[ \log(m_{ij}) = \lambda + \lambda_i^{JK} + \lambda_j^{Pref} + \lambda_{ij}^{JKxPref} \]

dengan:

$m_{ij}$ = nilai harapan frekuensi untuk sel ke-$(i,j)$
$\lambda$ = intercept (rata-rata log-count semua sel)
$\lambda_i^{\text{JK}}$ = efek baris (Jenis Kelamin)
$\lambda_j^{\text{Pref}}$ = efek kolom (Preferensi)
$\lambda_{ij}^{\text{JK} \times \text{Pref}}$ = efek interaksi

15.11 Fitting Model Log-Linear di R

# Membuat data frame dari tabel
 tabel2x3 <- matrix(c(15, 18, 7,
                        22, 20, 10), nrow = 2, byrow = TRUE)
 colnames(tabel2x3) <- c("Tatap Muka", "Daring", "Campuran")
 rownames(tabel2x3) <- c("Laki-laki", "Perempuan")
 tabel2x3

##           Tatap Muka Daring Campuran
## Laki-laki         15     18        7
## Perempuan         22     20       10

# Ubah menjadi data.frame untuk glm
 data2x3 <- as.data.frame(as.table(tabel2x3))
 colnames(data2x3) <- c("JenisKelamin", "Metode", "Freq")
 data2x3

##   JenisKelamin     Metode Freq
## 1    Laki-laki Tatap Muka   15
## 2    Perempuan Tatap Muka   22
## 3    Laki-laki     Daring   18
## 4    Perempuan     Daring   20
## 5    Laki-laki   Campuran    7
## 6    Perempuan   Campuran   10

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

## 
## Call:
## glm(formula = Freq ~ JenisKelamin + Metode, family = poisson, 
##     data = data2x3)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            2.77801    0.20287  13.693  < 2e-16 ***
## JenisKelaminPerempuan  0.26236    0.21031   1.248  0.21221    
## MetodeDaring           0.02667    0.23096   0.115  0.90808    
## MetodeCampuran        -0.77770    0.29300  -2.654  0.00795 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 12.09914  on 5  degrees of freedom
## Residual deviance:  0.40016  on 2  degrees of freedom
## AIC: 35.438
## 
## Number of Fisher Scoring iterations: 3

15.12 Intepretasi

Interpretasi Koefisien:

Koefisien	Estimasi	Arti
(Intercept)	2.778	Log-frekuensi untuk referensi: Laki-laki dengan Metode Tatap Muka.
JenisKelaminPerempuan	0.262	Perempuan cenderung memiliki frekuensi sekitar 30% lebih tinggi dari laki-laki karena tidak signifikan (p = 0.21).
MetodeDaring	0.0267	Frekuensi metode daring hampir sama dengan tatap muka (baseline), tidak signifikan (p = 0.91).
MetodeCampuran	-0.7777	Frekuensi metode campuran sekitar 54% lebih rendah dari tatap muka signifikan (p = 0.00795).

Goodness of Fit:

Residual Deviance = 0.40 dengan 2 degrees of freedom: Artinya model fit dengan baik ke data (karena deviance kecil).
AIC = 35.438: Digunakan untuk membandingkan antar model (misalnya, jika kamu membuat model dengan interaksi).

Kesimpulan

Model tanpa interaksi ini cocok digunakan (deviance kecil).

Tidak ada bukti yang signifikan bahwa jenis kelamin dan metode daring berpengaruh secara individual terhadap frekuensi.

Metode Campuran berbeda secara signifikan dari Tatap Muka dalam hal frekuensi responden.

Karena model ini tidak mengandung interaksi, kita mengasumsikan Jenis Kelamin dan Metode bersifat independen.

16 Model Log Linear Tiga Arah

Model log-linear tiga arah digunakan untuk menganalisis tabel kontingensi tiga dimensi, misalnya melibatkan tiga variabel kategorik: $A$, $B$, dan $C$

16.1 Model Log-Linear untuk Tabel Tiga Arah

1. Model Saturated (Penuh)

Model ini memuat semua efek utama, interaksi dua arah, dan tiga arah.

\[ \log(m_{ijk}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{ij}^{AB} + \lambda_{ik}^{AC} + \lambda_{jk}^{BC} + \lambda_{ijk}^{ABC} \]

2. Model Homogen Association

Model ini mengasumsikan bahwa semua interaksi dua arah signifikan, tetapi tidak ada interaksi tiga arah.

\[ \log(m_{ijk}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{ij}^{AB} + \lambda_{ik}^{AC} + \lambda_{jk}^{BC} \]

3. Model Conditional Independence

Misalnya, $A\perp B \mid$ C (A dan B independen kondisional terhadap C):

\[ \log(m_{ijk}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{ik}^{AC} + \lambda_{jk}^{BC} \]

4. Model Joint Independence

Misalnya, $A\perp (B, C)$— A independen dari gabungan B dan C:

\[ \log(m_{ijk}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{jk}^{BC} \]

5. Model Tanpa Interaksi (Mutual Independence)

Asumsi bahwa semua variabel saling independen (tidak ada interaksi sama sekali):

\[ \log(m_{ijk}) = \lambda + \lambda_i^A + \lambda_j^B + \lambda_k^C \]

16.2 Pengujian Interaksi dalam Model Log-Linear Tiga Arah

Pengujian interaksi dalam model log-linear tiga arah bertujuan untuk mengetahui apakah terdapat hubungan (asosiasi) yang signifikan antara ketiga variabel kategorik dalam sebuah tabel kontingensi tiga arah. Proses pengujian ini dilakukan dengan membandingkan model-model log-linear yang memiliki dan tidak memiliki komponen interaksi tertentu.

16.3 Analisis Log-Linear untuk Tabel Tiga Arah

Data tentang preferensi jenis transportasi berdasarkan usia dan waktu tempuh ke kampus.

 library("epitools")
 library("DescTools")
 library("lawstat")

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

# Variabel faktor
transport <- factor(rep(c("Motor", "Mobil", "Sepeda"), each = 4))
usia <- factor(rep(c("Remaja", "Dewasa"), each = 2, times = 3))
waktu <- factor(rep(c("Cepat", "Lama"), times = 6))
counts <- c(150, 60, 90, 40, 80, 45, 30, 25, 50, 20, 70, 35)

# Membuat data frame
data <- data.frame(
  Transportasi = transport,
  Kelompok_Usia = usia,
  Waktu_Tempuh = waktu,
  Frekuensi = counts
)

data

##    Transportasi Kelompok_Usia Waktu_Tempuh Frekuensi
## 1         Motor        Remaja        Cepat       150
## 2         Motor        Remaja         Lama        60
## 3         Motor        Dewasa        Cepat        90
## 4         Motor        Dewasa         Lama        40
## 5         Mobil        Remaja        Cepat        80
## 6         Mobil        Remaja         Lama        45
## 7         Mobil        Dewasa        Cepat        30
## 8         Mobil        Dewasa         Lama        25
## 9        Sepeda        Remaja        Cepat        50
## 10       Sepeda        Remaja         Lama        20
## 11       Sepeda        Dewasa        Cepat        70
## 12       Sepeda        Dewasa         Lama        35

Membentuk Tabel Kontigensi 3 Arah

table3d<-xtabs(Frekuensi~ Transportasi+ Kelompok_Usia + Waktu_Tempuh, data= data)
 ftable(table3d)

##                            Waktu_Tempuh Cepat Lama
## Transportasi Kelompok_Usia                        
## Mobil        Dewasa                        30   25
##              Remaja                        80   45
## Motor        Dewasa                        90   40
##              Remaja                       150   60
## Sepeda       Dewasa                        70   35
##              Remaja                        50   20

17 UJI MODEL INTERAKSI TIGA ARAH (SATURATED VS HOMOGENOUS)

17.0.1 Penentuan Kategori Referensi

 ##=============================##
 #Penentuankategorireference
 ##=============================##
usia <- relevel(usia, ref = "Remaja")
waktu <-relevel(waktu,ref="Lama")
 transport <-relevel(transport,ref="Sepeda")

#Modelsaturated
 model_saturated<-glm(counts ~usia +waktu + transport+
 usia*waktu +usia*transport+ waktu*transport+
 usia*waktu*transport,
 family=poisson(link ="log"))
 summary(model_saturated)

## 
## Call:
## glm(formula = counts ~ usia + waktu + transport + usia * waktu + 
##     usia * transport + waktu * transport + usia * waktu * transport, 
##     family = poisson(link = "log"))
## 
## Coefficients:
##                                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                           2.996e+00  2.236e-01  13.397  < 2e-16 ***
## usiaDewasa                            5.596e-01  2.803e-01   1.996 0.045885 *  
## waktuCepat                            9.163e-01  2.646e-01   3.463 0.000534 ***
## transportMobil                        8.109e-01  2.687e-01   3.018 0.002549 ** 
## transportMotor                        1.099e+00  2.582e-01   4.255 2.09e-05 ***
## usiaDewasa:waktuCepat                -2.231e-01  3.359e-01  -0.664 0.506542    
## usiaDewasa:transportMobil            -1.147e+00  3.752e-01  -3.058 0.002229 ** 
## usiaDewasa:transportMotor            -9.651e-01  3.468e-01  -2.783 0.005383 ** 
## waktuCepat:transportMobil            -3.409e-01  3.236e-01  -1.054 0.292104    
## waktuCepat:transportMotor            -1.641e-15  3.055e-01   0.000 1.000000    
## usiaDewasa:waktuCepat:transportMobil -1.699e-01  4.700e-01  -0.361 0.717743    
## usiaDewasa:waktuCepat:transportMotor  1.178e-01  4.151e-01   0.284 0.776601    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  2.2101e+02  on 11  degrees of freedom
## Residual deviance: -3.6415e-14  on  0  degrees of freedom
## AIC: 92.915
## 
## Number of Fisher Scoring iterations: 3

 exp(model_saturated$coefficients)

##                          (Intercept)                           usiaDewasa 
##                           20.0000000                            1.7500000 
##                           waktuCepat                       transportMobil 
##                            2.5000000                            2.2500000 
##                       transportMotor                usiaDewasa:waktuCepat 
##                            3.0000000                            0.8000000 
##            usiaDewasa:transportMobil            usiaDewasa:transportMotor 
##                            0.3174603                            0.3809524 
##            waktuCepat:transportMobil            waktuCepat:transportMotor 
##                            0.7111111                            1.0000000 
## usiaDewasa:waktuCepat:transportMobil usiaDewasa:waktuCepat:transportMotor 
##                            0.8437500                            1.1250000

17.0.2 Ringkasan Model

Model yang digunakan adalah model log-linear jenuh yang mencakup semua efek utama, interaksi dua variabel, serta interaksi tiga variabel. Model ini bertujuan untuk memodelkan hubungan antara Kelompok Usia, Waktu Tempuh, dan Transportasi berdasarkan frekuensi responden.

17.0.3 Hasil Estimasi Koefisien

Coefficients:                                  Estimate Std. Error (Intercept)                      4.248495   0.119523 usia1M                          -0.356675   0.186263 waktu1fav                        0.461035   0.152626 transport1fund                   0.041964   0.167285 transport2mod                    0.405465   0.154303 usia1M:waktu1fav                 0.426268   0.228268 usia1M:transport1fund           -0.468049   0.282210 usia1M:transport2mod            -0.271934   0.249148 waktu1fav:transport1fund         0.060690   0.212423 waktu1fav:transport2mod          0.008969   0.196903 usia1M:waktu1fav:transport1fund  0.438301   0.336151 usia1M:waktu1fav:transport2mod   0.282383   0.301553
                                z value Pr(>|z|)    
(Intercept)                      35.545  < 2e-16 ***
usia1M                           -1.915  0.05551 .  
waktu1fav                         3.021  0.00252 ** 
transport1fund                    0.251  0.80193    
transport2mod                     2.628  0.00860 ** 
usia1M:waktu1fav                  1.867  0.06185 .  
usia1M:transport1fund            -1.659  0.09721 .  
usia1M:transport2mod             -1.091  0.27507    
waktu1fav:transport1fund          0.286  0.77511    
waktu1fav:transport2mod           0.046  0.96367    
usia1M:waktu1fav:transport1fund   1.304  0.19227    
usia1M:waktu1fav:transport2mod    0.936  0.34905

17.0.4 Interpretasi Koefisien

Efek Utama

(Intercept): 4.248
Ini adalah log(count) untuk kategori referensi: usia=Dewasa, waktu=Lama, transport=Sepeda.
→ $\exp(4.248)$ $\approx$ 70: nilai rata-rata (frekuensi) untuk kategori referensi.
usia1M: -0.357 (p = 0.055)
Remaja cenderung memiliki penurunan frekuensi sebesar ~30% dibanding dewasa
(karena $\exp(-0.357)$ $\approx$ 70 , marginally signifikan.
waktu1fav: 0.461 (p = 0.0025)
Waktu tempuh cepat cenderung menaikkan frekuensi sebesar 58% dibanding waktu lama. (karena exp(0.461)≈1.59), signifikan.

**transport1fund: 0.042 (p = 0.802)**\
Tidak signifikan; tidak ada bukti perbedaan dengan referensi `Sepeda`.

transport2mod: 0.405 (p = 0.009)
Transportasi Mobil cenderung meningkatkan frekuensi sebesar 50% dibanding Sepeda, signifikan.

Efek Interaksi Dua Arah

usia1M:waktu1fav: 0.426 (p = 0.062)
Kombinasi remaja + waktu cepat cenderung meningkatkan frekuensi (sekitar 53%) dibanding ekspektasi berdasarkan efek terpisah, namun marginal.
usia1M:transport1fund: -0.468 (p = 0.097) Interaksi remaja + motor menunjukkan kecenderungan penurunan frekuensi (sekitar -37%), marginal.
Interaksi lainnya tidak signifikan secara statistik (p > 0.1).

Efek Interaksi Tiga Arah

Tidak ada yang signifikan (p > 0.1). Artinya, tidak ada bukti kuat bahwa hubungan antara dua variabel tergantung pada kategori ketiga secara signifikan.

🔹 Kesimpulan Sederhana

Variabel yang paling berpengaruh signifikan terhadap jumlah frekuensi adalah:

-   `waktu` (lebih cepat → frekuensi naik)

-   `transport` (mobil → frekuensi naik)

Beberapa interaksi dua arah marginal signifikan, terutama antara usia dan waktu.
Efek tiga arah tidak signifikan, sehingga model bisa dipertimbangkan untuk disederhanakan.

17.1 Model Homogenous

Model homogenous adalah model statistik yang mengasumsikan bahwa interaksi antara variabel-variabel dalam tabel kontingensi bersifat seragam atau konstan di seluruh kategori. Dengan kata lain, model ini menganggap bahwa efek interaksi tidak berbeda-beda antar kelompok, sehingga pola hubungan antar variabel dapat digambarkan dengan satu set parameter yang sama untuk seluruh data. Model ini biasanya digunakan untuk menyederhanakan analisis ketika interaksi yang kompleks dianggap tidak signifikan.

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

# Homogenous Model
 model_homogenous <- glm(counts ~ usia + waktu + transport +
 usia*waktu + usia*transport + waktu*transport,
 family = poisson(link = "log"))
 summary(model_homogenous)

## 
## Call:
## glm(formula = counts ~ usia + waktu + transport + usia * waktu + 
##     usia * transport + waktu * transport, family = poisson(link = "log"))
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.00104    0.18262  16.434  < 2e-16 ***
## usiaDewasa                 0.55127    0.19412   2.840  0.00451 ** 
## waktuCepat                 0.90886    0.19371   4.692 2.71e-06 ***
## transportMobil             0.84227    0.21607   3.898 9.69e-05 ***
## transportMotor             1.06309    0.20152   5.275 1.33e-07 ***
## usiaDewasa:waktuCepat     -0.21115    0.16899  -1.249  0.21149    
## usiaDewasa:transportMobil -1.24525    0.22453  -5.546 2.92e-08 ***
## usiaDewasa:transportMotor -0.88283    0.19065  -4.631 3.64e-06 ***
## waktuCepat:transportMobil -0.39136    0.22946  -1.706  0.08809 .  
## waktuCepat:transportMotor  0.04949    0.20521   0.241  0.80942    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 221.0056  on 11  degrees of freedom
## Residual deviance:   0.4958  on  2  degrees of freedom
## AIC: 89.411
## 
## Number of Fisher Scoring iterations: 3

17.2 Uji Hipotesis: Apakah Ada Interaksi Tiga Arah? (Saturated vs Homogenous)

Pengujian ini menggunakan residual deviance dari kedua model (saturated dan homogenous).

17.2.1 Langkah-Langkah Pengujian

17.2.1.1 1. Hipotesis

H0: Tidak ada interaksi tiga arah (model homogenous sudah cukup)
H1: Ada interaksi tiga arah (model saturated diperlukan)

17.2.1.2 2. Hitung Selisih Deviance

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

## [1] 0.4958032

17.2.1.3 3. Hitung Derajat Bebas

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

## [1] 2

17.2.1.4 4. Chi-Square Tabel (𝛼 = 0.05)

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

## [1] 5.991465

17.2.1.5 5. Keputusan Uji

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

## [1] "Terima H0"

Kesimpulan

Setelah dilakukan uji dengan alpha = 5%, disimpulkan bahwa tidak ada interaksi 3 arah

18 UJI MODEL INTERAKSI DUA ARAH(HOMOGENOUS VSCONDITIONAL ON X)

18.1 Model Conditional on X

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

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

# Conditional Association on X
 model_conditional_X <- glm(counts ~ usia + waktu + transport +
 usia*waktu + usia*transport,
 family = poisson(link = "log"))
 summary(model_conditional_X)

## 
## Call:
## glm(formula = counts ~ usia + waktu + transport + usia * waktu + 
##     usia * transport, family = poisson(link = "log"))
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 3.0729     0.1408  21.829  < 2e-16 ***
## usiaDewasa                  0.5163     0.1895   2.725 0.006423 ** 
## waktuCepat                  0.8065     0.1076   7.497 6.52e-14 ***
## transportMobil              0.5798     0.1493   3.884 0.000103 ***
## transportMotor              1.0986     0.1380   7.960 1.72e-15 ***
## usiaDewasa:waktuCepat      -0.1646     0.1638  -1.005 0.314925    
## usiaDewasa:transportMobil  -1.2264     0.2236  -5.485 4.13e-08 ***
## usiaDewasa:transportMotor  -0.8850     0.1904  -4.648 3.36e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 221.0056  on 11  degrees of freedom
## Residual deviance:   5.8883  on  4  degrees of freedom
## AIC: 90.804
## 
## Number of Fisher Scoring iterations: 4

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

18.2.1 Hipotesis

H0 : (Tidak ada interaksi antara Waktu Tempuh dan Transport

H1:(ada interaksi antara Waktu Tempuh dan Transport

18.2.2 Tingkat Signifikansi

𝛼=5%

18.2.3 Statistik Uji

Deviance = Deviance model conditional on X−Deviance model homogenous =3.903−1.798 = 2.132 𝑑𝑏=𝑑𝑏 model conditional on X−𝑑𝑏model homogenous =4−2=2

18.2.4 Daerah Penolakan

• Tolak 𝐻0 jika ΔDeviance >𝜒2 0.05,2 = 5.991

18.2.5 Keputusan

• Karena 2.132 < 5.991, maka tidak tolak 𝐻0

18.2.6 Kesimpulan

• Dengan taraf nyata 5%, belum cukup bukti untuk menolak 𝐻0 atau dapat dikatakan bahwa tidak ada interaksi antara Waktu Tempuh dan Transport

# Deviance of Model
 Deviance.model <- model_conditional_X$deviance- model_homogenous$deviance # model_conditional_X: conditional on X, model_homogenous: homogenous
 Deviance.model

## [1] 5.39249

18.3 Pengujian Selisih Deviance (Conditional on X vs Homogenous)

 # Selisih deviance antar model
 Deviance.model <- model_conditional_X$deviance- model_homogenous$deviance
 Deviance.model

## [1] 5.39249

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

## [1] 2

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

## [1] 5.991465

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

## [1] "Terima"

19 UJI MODEL INTERAKSI DUA ARAH(HOMOGENOUS VS CONDITIONAL ON Y)

19.1 Model Conditional on Y

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

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

# Conditional Association on Y
 model_conditional_Y <- glm(counts ~ usia + waktu + transport +
 usia*waktu + waktu*transport,
 family = poisson(link = "log"))
 summary(model_conditional_Y)

## 
## Call:
## glm(formula = counts ~ usia + waktu + transport + usia * waktu + 
##     waktu * transport, family = poisson(link = "log"))
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.41955    0.14744  23.193  < 2e-16 ***
## usiaDewasa                -0.22314    0.13416  -1.663 0.096269 .  
## waktuCepat                 0.85000    0.17752   4.788 1.68e-06 ***
## transportMobil             0.24116    0.18019   1.338 0.180767    
## transportMotor             0.59784    0.16787   3.561 0.000369 ***
## usiaDewasa:waktuCepat     -0.16462    0.16381  -1.005 0.314925    
## waktuCepat:transportMobil -0.32817    0.22336  -1.469 0.141771    
## waktuCepat:transportMotor  0.09531    0.20170   0.473 0.636541    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 221.006  on 11  degrees of freedom
## Residual deviance:  35.933  on  4  degrees of freedom
## AIC: 120.85
## 
## Number of Fisher Scoring iterations: 4

19.2 Pengujian Ada Tidaknya Interaksi antara X dan Z (Homogenous model vs Conditional Association on Y)

Hipotesis

– 𝐻0 ∶ (Tidak ada interaksi antara usia dan transport

– 𝐻1 ∶(Ada interaksi antara jusia dan transport
Tingkat Signifikansi

– 𝛼=5%
Statistik Uji

– ΔDeviance = Deviance model conditional on Y−Deviance model homogenous– =2.9203−1.798 = 1.1223

– 𝑑𝑏=𝑑𝑏 model conditional on Y−𝑑𝑏model homogenous– =4−2=2
Daerah Penolakan – Tolak 𝐻0 jika ΔDeviance >𝜒2 0.05,2 = 5.991
Keputusan– Karena 1.1223 < 5.991, maka tidak tolak 𝐻0
Kesimpulan– Dengan taraf nyata 5%, belum cukup bukti untuk menolak

19.3 Pengujian Hipotesis Interaksi X dan Z (Conditional on Y vs Homogenous)

# Deviance of Model
 Deviance.model <- model_conditional_Y$deviance- model_homogenous$deviance # model_conditional_Y: conditional on Y, model_homogenous: homogenous
 Deviance.model

## [1] 35.43675

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

## [1] 2

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

## [1] 5.991465

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

## [1] "Tolak"

Kesimpulan:

Pada taraf nyata 5%, belum cukup bukti untuk menolak 𝐻0. Artinya, tidak ada interaksi antara usia dan transport

20 UJI MODEL INTERAKSI DUA ARAH(HOMOGENOUS VS CONDITIONAL ON Z)

20.1 Model Conditional on Z

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

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

# Conditional Association on Z
 model_conditional_Z <- glm(counts ~ usia + waktu + transport +
 usia*transport + waktu*transport,
 family = poisson(link = "log"))
 summary(model_conditional_Z)

## 
## Call:
## glm(formula = counts ~ usia + waktu + transport + usia * transport + 
##     waktu * transport, family = poisson(link = "log"))
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.09104    0.16356  18.898  < 2e-16 ***
## usiaDewasa                 0.40547    0.15430   2.628 0.008596 ** 
## waktuCepat                 0.78016    0.16283   4.791 1.66e-06 ***
## transportMobil             0.79281    0.20853   3.802 0.000144 ***
## transportMotor             1.03229    0.19640   5.256 1.47e-07 ***
## usiaDewasa:transportMobil -1.22645    0.22359  -5.485 4.13e-08 ***
## usiaDewasa:transportMotor -0.88504    0.19043  -4.648 3.36e-06 ***
## waktuCepat:transportMobil -0.32817    0.22336  -1.469 0.141771    
## waktuCepat:transportMotor  0.09531    0.20170   0.473 0.636541    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 221.0056  on 11  degrees of freedom
## Residual deviance:   2.0548  on  3  degrees of freedom
## AIC: 88.97
## 
## Number of Fisher Scoring iterations: 4

20.2 Pengujian Ada Tidaknya Interaksi antara X dan Y (Homogenous model vs Conditional Association on Z)

Hipotesis

– 𝐻0 ∶ Tidak ada interaksi antara usia dan transport

– 𝐻1 ∶ Ada interaksi antara usia dan transport
Tingkat Signifikansi

– 𝛼=5%
Statistik Uji

– ΔDeviance = Deviance model conditional on Z

−Deviance model homogenous

– =29.729−1.798 = 27.931

– 𝑑𝑏=𝑑𝑏 model conditional on Z

−𝑑𝑏model homogenous– =3−2=1 Daerah Penolakan– Tolak 𝐻0 jika ΔDeviance >𝜒2 0.05,1 = 3.841
Keputusan– Karena 27.931 > 3.841, maka tolak 𝐻0
Kesimpulan

– Dengan taraf nyata 5 persen, ada interaksi antara jenis kelamin dan pendapat tentang hukuman mati.

20.3 Pengujian Hipotesis Interaksi X dan Y (Conditional on Z vs Homogenous)//

# Deviance of Model
 Deviance.model <- model_conditional_Z$deviance- model_homogenous$deviance # model_conditional_Z: conditional on Z, model_homogenous: homogenous
 Deviance.model

## [1] 1.559001

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

## [1] 1

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

## [1] 3.841459

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

## [1] "Terima"

Kesimpulan: Pada taraf nyata 5%, terdapat bukti yang cukup untuk menolak 𝐻0. Artinya, ada interaksi yang signifikan antara jenis kelamin (X) dan pendapat tentang hukuman mati(Y).

21 PEMILIHAN MODEL TERBAIK

21.1 Ringkasan Model Log Linier

Model	Jumlah Parameter	Deviance	Derajat Bebas (df)	AIC

Saturated

≈ 0

92.915

Homogeneous

0.496

89.411

Conditional on X

5.888

90.804

Conditional on Y

35.933

120.850

Conditional on Z

2.055

88.970

Kesimpulan

Model Saturated cocok sempurna (deviance ≈ 0) tapi terlalu kompleks (df = 0).
Model Homogeneous memberikan fit yang sangat baik dengan kompleksitas lebih rendah dari saturated.
Model Conditional on Z juga menunjukkan fit baik (deviance rendah, AIC terendah)
Model Conditional on Y buruk karena deviance tinggi dan AIC besar → fit buruk.
Model Conditional on X cukup baik, tapi tidak sebaik Homogeneous dan Z.

21.2 Ringkasan Pengujian Interaksi 3 Arah dan 2 Arah

Pengujian	Deviance	df	χ² tabel (α=0.05)	Keputusan	Keterangan
XYZ vs YZ,XZ,XY (uji interaksi 3 arah)	0.496	2	5.991	Gagal tolak H₀	Interaksi 3 arah tidak signifikan
YZ vs XZ,XY (uji interaksi YZ)	33.878	2	5.991	Tolak H₀	Interaksi YZ signifikan
XZ vs YZ,XY (uji interaksi XZ)	1.559	1	3.841	Gagal tolak H₀	Interaksi XZ tidak signifikan
XY vs YZ,XZ (uji interaksi XY)	3.833	1	3.841	Gagal/tolak H₀ (batas)	Interaksi XY mendekati signifikan (marginal)

21.3 Kesimpulan Pemilihan Model Terbaik

Dari hasil di atas diketahui bahwa asosiasi yang nyata hanya terdapat interaksi antara waktu dan transport Sehingga, model terbaik adalah:

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

Model terbaik adalah model log-linear tanpa interaksi tiga arah dan hanya memuat interaksi dua arah antara kelompok usia dan transport.

22 MODEL TERBAIK

Model terbaik dipilih berdasarkan pengujian interaksi yang signifikan, yaitu hanya interaksi dua arah antara waktu dan transport

# Model Terbaik
 bestmodel <- glm(counts ~ usia + waktu + transport +
 transport*waktu,
 family = poisson(link = "log"))
 summary(bestmodel)

## 
## Call:
## glm(formula = counts ~ usia + waktu + transport + transport * 
##     waktu, family = poisson(link = "log"))
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                3.46731    0.13861  25.015  < 2e-16 ***
## usiaDewasa                -0.33401    0.07692  -4.342 1.41e-05 ***
## waktuCepat                 0.78016    0.16283   4.791 1.66e-06 ***
## transportMobil             0.24116    0.18019   1.338 0.180766    
## transportMotor             0.59784    0.16787   3.561 0.000369 ***
## waktuCepat:transportMobil -0.32817    0.22336  -1.469 0.141770    
## waktuCepat:transportMotor  0.09531    0.20170   0.473 0.636540    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 221.01  on 11  degrees of freedom
## Residual deviance:  36.94  on  5  degrees of freedom
## AIC: 119.86
## 
## Number of Fisher Scoring iterations: 4

23 INTERPRETASI KOEFISIEN MODEL TERBAIK

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

##                                  koef   exp_koef
## (Intercept)                3.46730841 32.0503597
## usiaDewasa                -0.33400614  0.7160494
## waktuCepat                 0.78015856  2.1818182
## transportMobil             0.24116206  1.2727273
## transportMotor             0.59783700  1.8181818
## waktuCepat:transportMobil -0.32817343  0.7202381
## waktuCepat:transportMotor  0.09531018  1.1000000

23.1 Interpretasi Koefisien Model Terbaik

Interpretasi Tiap Variabel

Variabel	Estimate	Interpretasi (log)	Interpretasi (eksp)
(Intercept)	3.467	Log rata-rata counts untuk usia Anak, waktu Lambat, transportasi Sepeda.	exp(3.467) ≈ 32.05 → jumlah rata-rata 32 responden.
usiaDewasa	-0.334	Dewasa cenderung memiliki log-count lebih rendah 0.334 dari Anak.	exp(-0.334) ≈ 0.716 → jumlahnya sekitar 71.6% dari anak-anak.
waktuCepat	0.780	Waktu cepat menaikkan log-count sebesar 0.78 dibanding waktu lambat.	exp(0.780) ≈ 2.18 → jumlahnya 2.18 kali lipat lebih banyak.
transportMobil	0.241	Mobil (vs Sepeda) meningkatkan log-count, tapi tidak signifikan.	exp(0.241) ≈ 1.27 → 27% lebih banyak, tidak signifikan (p = 0.18).
transportMotor	0.598	Motor (vs Sepeda) secara signifikan meningkatkan log-count sebesar 0.598.	exp(0.598) ≈ 1.82 → jumlahnya 1.82 kali lipat, signifikan.
waktuCepat:transportMobil	-0.328	Ada penurunan tambahan pada log-count untuk kombinasi waktu cepat dan mobil.	exp(-0.328) ≈ 0.72 → efek interaksi menurunkan hasil sekitar 28%, tidak signifikan (p = 0.14).
waktuCepat:transportMotor	0.095	Kombinasi waktu cepat dan motor sedikit menaikkan log-count, tapi tidak signifikan.	exp(0.095) ≈ 1.10 → hanya 10% peningkatan, tidak signifikan (p = 0.64).

Kesimpulan Utama:

Usia Dewasa → Signifikan menurunkan counts (responden lebih sedikit).
Waktu Cepat → Sangat signifikan meningkatkan counts.
Transportasi Motor → Signifikan meningkatkan counts (dibanding sepeda).
Transportasi Mobil dan interaksinya → Tidak signifikan.
Interaksi Waktu × Transportasi → Tidak signifikan secara statistik, tapi dipertahankan karena model terbaik berdasarkan uji sebelumnya.

24 NILAI DUGAAN MODEL TERBAIK

# Fitted values dari model terbaik
 data.frame(
 Transportasi = transport,
 Kelompok_Usia = usia,
 Waktu_Tempuh = waktu,
 Frekuensi = counts,
 fitted = bestmodel$fitted.values
 )

##    Transportasi Kelompok_Usia Waktu_Tempuh Frekuensi    fitted
## 1         Motor        Remaja        Cepat       150 139.85612
## 2         Motor        Remaja         Lama        60  58.27338
## 3         Motor        Dewasa        Cepat        90 100.14388
## 4         Motor        Dewasa         Lama        40  41.72662
## 5         Mobil        Remaja        Cepat        80  64.10072
## 6         Mobil        Remaja         Lama        45  40.79137
## 7         Mobil        Dewasa        Cepat        30  45.89928
## 8         Mobil        Dewasa         Lama        25  29.20863
## 9        Sepeda        Remaja        Cepat        50  69.92806
## 10       Sepeda        Remaja         Lama        20  32.05036
## 11       Sepeda        Dewasa        Cepat        70  50.07194
## 12       Sepeda        Dewasa         Lama        35  22.94964

24.1 Perhitungan Manual Nilai Dugaan (Fitted Value) Model Terbaik

No	Kombinasi	Frekuensi	Fitted	Interpretasi
1	Motor, Remaja, Cepat	150	139.86	Model cukup akurat, hanya selisih ~10 responden.
2	Motor, Remaja, Lama	60	58.27	Sangat akurat, selisih < 2.
3	Motor, Dewasa, Cepat	90	100.14	Prediksi agak berlebih (~10 responden).
4	Motor, Dewasa, Lama	40	41.73	Akurat.
5	Mobil, Remaja, Cepat	80	64.10	Model meremehkan responden (~16 lebih sedikit).
6	Mobil, Remaja, Lama	45	40.79	Akurat.
7	Mobil, Dewasa, Cepat	30	45.90	Model berlebihan (~16 responden).
8	Mobil, Dewasa, Lama	25	29.21	Cukup akurat.
9	Sepeda, Remaja, Cepat	50	69.93	Model berlebihan ~20 responden.
10	Sepeda, Remaja, Lama	20	32.05	Model berlebihan, prediksi lebih tinggi dari observasi.

Model mampu menangkap pola umum dari data cukup baik.
Selisih terbesar terjadi pada:
- Mobil, Remaja, Cepat (selisih -16)
- Mobil, Dewasa, Cepat (selisih +16)
- Sepeda, Remaja, Cepat (selisih +20)
Namun secara keseluruhan, prediksi mendekati frekuensi aktual → menunjukkan model cukup layak

25. Referensi

Fisher, R. A. (1935). The Design of Experiments. Oliver and Boyd.
Agresti, A. (2013). Categorical Data Analysis. John Wiley & Sons.
Mehta, C. R., & Patel, N. R. (1983). A Network Algorithm for Performing Fisher’s Exact Test in r × c Contingency Tables. Journal of the American Statistical Association, 78(382), 427-434.
Howell, D. C. (2012). Statistical Methods for Psychology (8th ed.). Cengage Learning.
Anderson, C. J. (2018). Log‑linear Models for Contingency Tables
Fox, J. (2008). Applied Regression Analysis and Generalized Linear Models.
Christensen, R. (1997). Log-Linear Models and Logistic Regression. Springer.
Bishop, Y. M. M., Fienberg, S. E., & Holland, P. W. (1975). Discrete Multivariate Analysis: Theory and Practice. MIT Press.
Holt, D. (1979). “Log‑Linear Models for Contingency Table Analysis.” Sociological Methods & Research.
Categorical Data Analysis: Chen, D., & Anderson, C. (2024). “Categorical Data Analysis”
Khamis, H. J. (1983). “Log‑linear model analysis of the semi‑symmetric intraclass contingency table.” Communications in Statistics – Theory and Methods.
Benoit, K. (2011). Multinomial and Ordinal Logistic Regression
Lee, E. (2018). Ordinal Logistic Regression and its Assumptions — Brant Test.

Cell	Keterangan	Nilai	Log(Frek)
\(\mu_{11}\)	Merokok = Ya, Sakit	30	log(30)≈3.401
\(\mu_{12}\)	Merokok = Ya, Sehat	20	log(20)≈2.996
\(\mu_{21}\)	Merokok = Tidak, Sakit	10	log(10)≈2.303
\(\mu_{22}\)	Merokok = Tidak, Sehat	40	log(40)≈3.689

E Book

Marhan Aristoteles Simanjuntak

2025-04-08

1. Analisis Data Kategori

1.1 Pendahuluan

1.2 Tujuan

2. Metode dalam Analisis Data Kategori

2.1 Tabel Kontingensi dan Uji Chi Square

2.2 Regresi Logistik (Logistic Regression)

3. Distribusi Probabilitas dalam Data Kategori

3.1 Distribusi Bernoulli

3.2 Distribusi Multinomial

3.3 Distribusi Binomial

3.4. Distribusi Poisson

4. Desain Sampling dalam analisis data kategori

4.1. Prospective Sampling

4.1.1. Studi Eksperimental

4.1.2. Studi Kohort

4.2. Retrospective Sampling

4.2.1. Studi Kasus-Kontrol

4.2.2. Studi Kohort Retrospektif

4.3 Tabel Perbandingan

5. Tabel Kontigensi 2x2

5.1 Distribusi Peluang dalam Tabel Kontigensi

5.1.1 Peluang Bersama (Joint Probability)

5.1.2 Peluang Marginal (Marginal Probability)

5.1.3 Peluang Bersyarat (Conditional Probability)

5.2 Ukuran Asosiasi dalam Data Kategori 2 x 2

5.2.1 Risk Difference (RD)

5.2.2 Relative Risk (RR)

5.2.2 Odds Ratio (OR)

6. Inferensi Tabel Kontigensi Dua Arah

6.1 Estimasi

6.1.1 Estimasi Titik

6.1.2 Estimasi Interval

6.2 Uji Hipotesis

6.2.1 Uji Proporsi

6.2.2 Uji Asosiasi

6.2.3 Uji Independensi

6.2.3.1 Uji Chi-Square

6.2.3.2 Uji Likelihood Ratio

6.2.3.3 Uji Exact Fisher

6.3 Analisis Residual dalam Tabel Kontigensi

7. Tabel Kontigensi Tiga Arah

7.1 Tabel Parsial dan Marginal

7.2 Distribusi Peluang

7.3 Tabel Peluang Bersyarat

7.4 Ukuran Asosiasi

7.4.1 Tabel Kontigensi Parsial

7.5 Conditional Independence

7.6 Inferensi Tabel Kontigensi Tiga Arah

7.6.1 Independensi Bersyarat dalam Tabel Kontigensi Tiga Arah

7.6.2 Pengujian Statistik untuk Independensi Bersyarat

7.6.3 Odds Ratio Bersama

7.6.4 Uji Homogenitas Odds Ratio dengan Statistik Breslow Day

8. Generalized Linear Model (GLM)

8.1 Exponential Family

8.2 Model Regresi Logistik

8.3 Model Regresi Poisson

9. Inferensi GLM

9.1 Mencari Ekspetasi dan Varians dalam GLM

9.2 Metode Penaksiran Parameter

9.3 Diagnostik Model GLM

9.4 Detail Metode Estimasi dan Inferensi Regresi Logistik

9.5 Detail Metode Estimasi dan Inferensi Regresi Poisson

10. Regresi Logistik dengan Prediktor Normal, Ordinal, dan Rasio

10.1 Simulasi Data

10.2 Eksplorasi Data

10.3 Perlakuan Variabel Ordinal

10.3.1 Treat sebagai Nominal (Dummy)

10.3.2 Treat sebagai Rasio

11. Pemilihan Model Regresi Logistik dan Evaluasi

11.1 Membangun Model Regresi Logistik : Pendekatan Confirmatory dan Exploratory

11.2 Metode Stepwise: Forward, Backward, dan Kedua Arah

11.3 Evaluasi Model ROC dan AUC

11.4 Pseudo R-Squared

11.5 Tabel Klasifikasi dan Evaluasi

11.6 Metode Perbandingan Model dalam Regresi Logistik

11.7 Likelihood-Ratio Test

11.8 Prinsip Parsimony