1 Pendahuluan

Regresi adalah salah satu alat inferensi statistik yang paling banyak digunakan untuk memodelkan hubungan antara variabel respons dan satu atau lebih variabel prediktor. Ketika variabel respons bukan berupa data kontinu, dibutuhkan pendekatan regresi khusus yang disesuaikan dengan skala pengukurannya.

Laporan ini membahas empat jenis regresi untuk data kategorik dan cacahan:

Jenis Regresi	Variabel Respons	Link Function	Package	Dataset
Biner	2 kategori (0/1)	Logit	`base R`	`Pima.tr` (MASS)
Multinomial	≥ 3 kategori tak terurut	Logit generalisasi	`nnet`	`hsbdemo` (UCLA)
Ordinal	≥ 3 kategori terurut	Cumulative Logit	`ordinal`	`wine` (ordinal)
Poisson	Data cacahan (count)	Log	`base R` + `AER`	`Affairs` (AER)

Prinsip umum: Semua model di atas adalah bagian dari Generalized Linear Model (GLM), yang menghubungkan nilai harapan \(E(Y)\) dengan prediktor melalui sebuah link function \(g(\cdot)\), sehingga \(g(E(Y)) = \mathbf{X}\boldsymbol{\beta}\).

2 Instalasi dan Pemuatan Package

pkgs <- c("nnet", "MASS", "ordinal", "AER", "ggplot2",
          "dplyr", "tidyr", "knitr", "effects")
new_pkgs <- pkgs[!(pkgs %in% installed.packages()[, "Package"])]
if (length(new_pkgs)) install.packages(new_pkgs, quiet = TRUE)

library(nnet);    library(MASS);    library(ordinal)
library(AER);     library(ggplot2); library(dplyr)
library(tidyr);   library(knitr);   library(effects)

3 Regresi Logistik Biner

3.1 Konsep dan Landasan Teori

Regresi logistik biner digunakan ketika variabel respons bersifat dikotom, yaitu hanya memiliki dua kemungkinan nilai: sukses (1) atau gagal (0). Karena \(P(Y=1)\) adalah probabilitas yang berada di interval \([0,1]\), kita tidak bisa langsung menggunakan regresi linear biasa.

Solusinya adalah mentransformasi probabilitas tersebut menggunakan fungsi logit:

\[\text{logit}(p) = \log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p \tag{1}\]

sehingga probabilitas prediksi diperoleh dari fungsi sigmoid:

\[\hat{p} = P(Y = 1 \mid \mathbf{x}) = \frac{e^{\beta_0 + \beta_1 x_1 + \cdots}}{1 + e^{\beta_0 + \beta_1 x_1 + \cdots}} \tag{2}\]

3.1.1 Notasi

Simbol	Keterangan
\(p\)	Probabilitas kejadian (\(Y = 1\))
\(\text{odds} = p/(1-p)\)	Rasio peluang terjadi vs tidak terjadi
\(\beta_k\)	Koefisien log-odds untuk prediktor \(x_k\)
\(e^{\beta_k}\)	Odds Ratio (OR) — perubahan odds per satu satuan \(x_k\)

3.1.2 Interpretasi Odds Ratio

\[OR_k = e^{\beta_k} \tag{3}\]

\(OR = 1\) : prediktor tidak berpengaruh
\(OR > 1\) : meningkatkan odds kejadian
\(OR < 1\) : menurunkan odds kejadian

3.2 Dataset: `Pima.tr` — Diabetes Suku Pima

Dataset berisi data kesehatan 200 wanita suku Indian Pima, dengan variabel respons type (Yes = menderita diabetes, No = tidak). Variabel prediktor yang digunakan: glu (kadar glukosa plasma), bmi (indeks massa tubuh), age (usia), dan bp (tekanan darah diastolik).

data(Pima.tr, package = "MASS")
kable(head(Pima.tr, 8), caption = "8 Baris Pertama Dataset Pima.tr")

8 Baris Pertama Dataset Pima.tr
npreg	glu	bp	skin	bmi	ped	age	type
5	86	68	28	30.2	0.364	24	No
7	195	70	33	25.1	0.163	55	Yes
5	77	82	41	35.8	0.156	35	No
0	165	76	43	47.9	0.259	26	No
0	107	60	25	26.4	0.133	23	No
5	97	76	27	35.6	0.378	52	Yes
3	83	58	31	34.3	0.336	25	No
1	193	50	16	25.9	0.655	24	No

cat("Dimensi data:", nrow(Pima.tr), "baris x", ncol(Pima.tr), "kolom\n")

## Dimensi data: 200 baris x 8 kolom

cat("\nDistribusi variabel respons (type):\n")

## 
## Distribusi variabel respons (type):

print(table(Pima.tr$type))

## 
##  No Yes 
## 132  68

cat("\nProporsi (%):\n")

## 
## Proporsi (%):

print(round(prop.table(table(Pima.tr$type)) * 100, 1))

## 
##  No Yes 
##  66  34

Catatan: Dataset ini memiliki distribusi kelas yang tidak seimbang. Proporsi kelas minoritas (Yes = diabetes) perlu diperhatikan saat mengevaluasi performa model.

3.3 Pemodelan

model_bin <- glm(type ~ glu + bmi + age + bp,
                 family = binomial(link = "logit"),
                 data   = Pima.tr)
summary(model_bin)

## 
## Call:
## glm(formula = type ~ glu + bmi + age + bp, family = binomial(link = "logit"), 
##     data = Pima.tr)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -9.144169   1.647351  -5.551 2.84e-08 ***
## glu          0.031224   0.006560   4.760 1.94e-06 ***
## bmi          0.093564   0.032645   2.866  0.00416 ** 
## age          0.054793   0.018166   3.016  0.00256 ** 
## bp          -0.006076   0.017345  -0.350  0.72612    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 256.41  on 199  degrees of freedom
## Residual deviance: 188.27  on 195  degrees of freedom
## AIC: 198.27
## 
## Number of Fisher Scoring iterations: 5

3.4 Odds Ratio dan Confidence Interval

or_table <- data.frame(
  Variabel   = names(coef(model_bin)),
  Koefisien  = round(coef(model_bin), 4),
  Odds_Ratio = round(exp(coef(model_bin)), 4),
  CI_Lower   = round(exp(confint(model_bin))[, 1], 4),
  CI_Upper   = round(exp(confint(model_bin))[, 2], 4)
)
kable(or_table, caption = "Odds Ratio dan Confidence Interval 95% — Model Biner")

Odds Ratio dan Confidence Interval 95% — Model Biner
	Variabel	Koefisien	Odds_Ratio	CI_Lower	CI_Upper
(Intercept)	(Intercept)	-9.1442	0.0001	0.0000	0.0023
glu	glu	0.0312	1.0317	1.0191	1.0458
bmi	bmi	0.0936	1.0981	1.0317	1.1734
age	age	0.0548	1.0563	1.0201	1.0959
bp	bp	-0.0061	0.9939	0.9603	1.0285

Interpretasi:

glu: Setiap kenaikan 1 unit kadar glukosa meningkatkan odds diabetes sebesar faktor \(e^{\hat\beta}\), dengan asumsi variabel lain tetap. Variabel ini konsisten menjadi prediktor terkuat diabetes.
bmi: Indeks massa tubuh yang lebih tinggi juga meningkatkan odds diabetes — sesuai fakta klinis bahwa obesitas merupakan faktor risiko utama.
age: Pertambahan usia meningkatkan odds diabetes secara signifikan.
bp: Tekanan darah memberikan pengaruh lebih kecil — perhatikan apakah \(p\)-value > 0.05 dan interval kepercayaan mencakup angka 1.

3.5 Visualisasi: Kurva Logistik

ggplot(Pima.tr, aes(x = glu,
                    y = as.numeric(type == "Yes"),
                    color = type)) +
  geom_point(alpha = 0.4, size = 2) +
  stat_smooth(method = "glm",
              method.args = list(family = "binomial"),
              se = TRUE, color = "#2c3e50", linewidth = 1.1) +
  scale_color_manual(values = c("No" = "#27ae60", "Yes" = "#c0392b")) +
  labs(
    title    = "Regresi Logistik Biner: Diabetes vs Kadar Glukosa",
    subtitle = "Kurva sigmoid menunjukkan probabilitas diabetes meningkat seiring glukosa",
    x     = "Kadar Glukosa / glu (mg/dL)",
    y     = "Probabilitas Diabetes P(Y=1)",
    color = "Status Diabetes"
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

3.6 Evaluasi Model

prob_pred  <- predict(model_bin, type = "response")
pred_class <- ifelse(prob_pred > 0.5, "Yes", "No")

cm <- table(Aktual = Pima.tr$type, Prediksi = pred_class)
kable(cm, caption = "Confusion Matrix — Regresi Biner")

Confusion Matrix — Regresi Biner
	No	Yes
No	113	19
Yes	30	38

akurasi      <- mean(pred_class == Pima.tr$type)
sensitivitas <- cm["Yes", "Yes"] / sum(cm["Yes", ])
spesifisitas <- cm["No",  "No"]  / sum(cm["No",  ])

cat("\nAkurasi     :", round(akurasi     * 100, 2), "%\n")

## 
## Akurasi     : 75.5 %

cat("Sensitivitas:", round(sensitivitas * 100, 2), "%\n")

## Sensitivitas: 55.88 %

cat("Spesifisitas:", round(spesifisitas * 100, 2), "%\n")

## Spesifisitas: 85.61 %

Interpretasi evaluasi: Akurasi saja tidak cukup untuk data tidak seimbang. Sensitivitas (kemampuan model mendeteksi kasus positif/diabetes) dan spesifisitas (kemampuan mendeteksi negatif) sama pentingnya, terutama dalam konteks medis di mana false negative berisiko tinggi.

4 Regresi Multinomial

4.1 Konsep dan Landasan Teori

Regresi multinomial merupakan generalisasi regresi logistik biner untuk situasi di mana variabel respons memiliki lebih dari dua kategori yang tidak memiliki urutan (nominal polytomous). Satu kategori dipilih sebagai referensi, dan model membandingkan setiap kategori lainnya terhadap referensi tersebut.

Untuk \(K\) kategori dengan kategori ke-\(K\) sebagai referensi:

\[\log\left(\frac{P(Y = k)}{P(Y = K)}\right) = \beta_{0k} + \beta_{1k} x_1 + \cdots + \beta_{pk} x_p, \quad k = 1, \ldots, K-1 \tag{4}\]

Probabilitas untuk setiap kategori diperoleh melalui fungsi softmax:

\[P(Y = k \mid \mathbf{x}) = \frac{e^{\eta_k}}{1 + \sum_{j=1}^{K-1} e^{\eta_j}} \tag{5}\]

4.1.1 Notasi

Simbol	Keterangan
\(K\)	Jumlah kategori respons
\(\eta_k = \mathbf{x}^T \boldsymbol{\beta}_k\)	Skor linear untuk kategori ke-\(k\)
\(e^{\beta_{jk}}\)	OR kategori ke-\(k\) vs referensi untuk prediktor \(x_j\)

4.2 Dataset: `hsbdemo` — Program Pendidikan Siswa

Dataset berisi data 200 siswa sekolah menengah dengan variabel respons prog: program pendidikan yang diikuti (academic, general, vocation). Kategori referensi yang dipilih adalah academic.

url <- "https://stats.idre.ucla.edu/stat/data/hsbdemo.csv"
hsbdemo <- tryCatch(
  read.csv(url),
  error = function(e) {
    set.seed(42)
    n <- 200
    data.frame(
      id      = 1:n,
      female  = sample(c(0, 1), n, replace = TRUE),
      ses     = sample(c("low", "middle", "high"), n, replace = TRUE),
      schtyp  = sample(c(1, 2), n, replace = TRUE),
      prog    = sample(c("general", "academic", "vocation"), n,
                       replace = TRUE, prob = c(0.2, 0.5, 0.3)),
      read    = round(rnorm(n, 52, 10)),
      write   = round(rnorm(n, 52, 10)),
      math    = round(rnorm(n, 52,  9)),
      science = round(rnorm(n, 51,  9)),
      socst   = round(rnorm(n, 52, 10)),
      honors  = sample(c("enrolled", "not enrolled"), n, replace = TRUE),
      awards  = sample(0:3, n, replace = TRUE),
      cid     = sample(1:20, n, replace = TRUE)
    )
  }
)

hsbdemo$prog <- relevel(as.factor(hsbdemo$prog), ref = "academic")

kable(head(hsbdemo[, c("prog", "ses", "write", "read", "math")], 8),
      caption = "8 Baris Pertama Dataset hsbdemo (variabel terpilih)")

8 Baris Pertama Dataset hsbdemo (variabel terpilih)
prog	ses	write	read	math
vocation	low	35	34	41
general	middle	33	34	41
vocation	high	39	39	44
vocation	low	37	37	42
vocation	middle	31	39	40
general	high	36	42	42
vocation	middle	36	31	46
vocation	middle	31	50	40

cat("Distribusi variabel respons (prog):\n")

## Distribusi variabel respons (prog):

print(table(hsbdemo$prog))

## 
## academic  general vocation 
##      105       45       50

cat("\nProporsi (%):\n")

## 
## Proporsi (%):

print(round(prop.table(table(hsbdemo$prog)) * 100, 1))

## 
## academic  general vocation 
##     52.5     22.5     25.0

4.3 Pemodelan

model_mn <- multinom(prog ~ ses + write + read, data = hsbdemo)

## # weights:  18 (10 variable)
## initial  value 219.722458 
## iter  10 value 176.897700
## final  value 175.171786 
## converged

summary(model_mn)

## Call:
## multinom(formula = prog ~ ses + write + read, data = hsbdemo)
## 
## Coefficients:
##          (Intercept)    seslow sesmiddle       write        read
## general     2.768623 0.9748410 0.5583957 -0.03179768 -0.04552660
## vocation    6.165563 0.6971476 1.1818269 -0.07451237 -0.07573009
## 
## Std. Errors:
##          (Intercept)    seslow sesmiddle      write       read
## general     1.375250 0.5262627 0.4709355 0.02494601 0.02398502
## vocation    1.441722 0.6119653 0.5226723 0.02542673 0.02666881
## 
## Residual Deviance: 350.3436 
## AIC: 370.3436

4.4 Koefisien, Z-Score, dan Odds Ratio

z_score <- summary(model_mn)$coefficients /
           summary(model_mn)$standard.errors
p_value <- (1 - pnorm(abs(z_score), 0, 1)) * 2

cat("=== Z-Score ===\n")

## === Z-Score ===

print(round(z_score, 3))

##          (Intercept) seslow sesmiddle  write   read
## general        2.013  1.852     1.186 -1.275 -1.898
## vocation       4.277  1.139     2.261 -2.930 -2.840

cat("\n=== P-Value ===\n")

## 
## === P-Value ===

print(round(p_value, 4))

##          (Intercept) seslow sesmiddle  write   read
## general       0.0441 0.0640    0.2357 0.2024 0.0577
## vocation      0.0000 0.2546    0.0238 0.0034 0.0045

cat("\n=== Odds Ratio ===\n")

## 
## === Odds Ratio ===

print(round(exp(coef(model_mn)), 3))

##          (Intercept) seslow sesmiddle write  read
## general       15.937  2.651     1.748 0.969 0.955
## vocation     476.069  2.008     3.260 0.928 0.927

Interpretasi:

Koefisien pada baris general dan vocation masing-masing dibandingkan terhadap kategori referensi academic.
Jika OR untuk write pada baris vocation < 1, artinya siswa dengan skor menulis lebih tinggi lebih kecil kemungkinannya memilih program vokasi dibandingkan program akademik.
Variabel yang memiliki \(|z| > 1.96\) dianggap signifikan secara statistik pada \(\alpha = 0.05\).
Fungsi multinom() dari package nnet tidak menampilkan p-value secara langsung, sehingga perlu dihitung manual dari z-score.

4.5 Visualisasi Prediksi

new_data_mn <- expand.grid(
  write = seq(min(hsbdemo$write), max(hsbdemo$write), length.out = 50),
  read  = mean(hsbdemo$read),
  ses   = "middle"
)
pred_prob <- predict(model_mn, newdata = new_data_mn, type = "probs")
plot_df   <- cbind(new_data_mn, as.data.frame(pred_prob)) %>%
  pivot_longer(cols     = colnames(as.data.frame(pred_prob)),
               names_to = "program", values_to = "probabilitas")

ggplot(plot_df, aes(x = write, y = probabilitas, color = program)) +
  geom_line(linewidth = 1.3) +
  scale_color_manual(values = c("academic" = "#2980b9",
                                "general"  = "#c0392b",
                                "vocation" = "#27ae60")) +
  labs(
    title    = "Probabilitas Prediksi Regresi Multinomial",
    subtitle = "Berdasarkan skor menulis — ses = middle, read = rata-rata",
    x     = "Skor Menulis (write)",
    y     = "Probabilitas",
    color = "Program"
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

Membaca grafik: Titik di mana kurva academic mulai menurun dan kurva vocation naik menunjukkan skor menulis ambang batas di mana program yang diprediksi berubah. Semakin tinggi skor menulis, semakin besar kecenderungan ke program academic.

5 Regresi Ordinal

5.1 Konsep dan Landasan Teori

Regresi ordinal (Cumulative Link Model / Proportional Odds Model) digunakan ketika variabel respons memiliki kategori yang terurut — misalnya penilaian kepuasan (rendah < sedang < tinggi). Model ini memanfaatkan informasi urutan yang tidak bisa ditangkap oleh regresi multinomial biasa.

Model Proportional Odds:

\[\log\left(\frac{P(Y \leq j)}{P(Y > j)}\right) = \alpha_j - (\beta_1 x_1 + \cdots + \beta_p x_p), \quad j = 1, \ldots, J-1 \tag{6}\]

di mana \(\alpha_j\) adalah threshold (intercept kumulatif) untuk setiap titik batas kategori.

5.1.1 Asumsi Proportional Odds

Asumsi kunci: Efek prediktor \(\boldsymbol{\beta}\) diasumsikan sama di semua titik batas kategori (parallel regression assumption). Asumsi ini wajib diuji menggunakan likelihood ratio test.

5.1.2 Notasi

Simbol	Keterangan
\(\alpha_j\)	Threshold (intercept) ke-\(j\)
\(P(Y \leq j)\)	Probabilitas kumulatif sampai kategori ke-\(j\)
\(e^{\beta_k}\)	Cumulative Odds Ratio untuk prediktor \(x_k\)

5.2 Dataset: `wine` — Penilaian Kepahitan Anggur

Dataset berisi 72 penilaian anggur oleh panel ahli. Variabel respons adalah rating kepahitan pada skala ordinal 1–5. Prediktor: temp (suhu penyimpanan: cold/warm) dan contact (kontak kulit: no/yes).

data(wine, package = "ordinal")
kable(head(wine, 8), caption = "8 Baris Pertama Dataset wine")

8 Baris Pertama Dataset wine
response	rating	temp	contact	bottle	judge
36	2	cold	no	1	1
48	3	cold	no	2	1
47	3	cold	yes	3	1
67	4	cold	yes	4	1
77	4	warm	no	5	1
60	4	warm	no	6	1
83	5	warm	yes	7	1
90	5	warm	yes	8	1

cat("Distribusi rating kepahitan (1 = paling tidak pahit, 5 = paling pahit):\n")

## Distribusi rating kepahitan (1 = paling tidak pahit, 5 = paling pahit):

print(table(wine$rating))

## 
##  1  2  3  4  5 
##  5 22 26 12  7

cat("\nProporsi (%):\n")

## 
## Proporsi (%):

print(round(prop.table(table(wine$rating)) * 100, 1))

## 
##    1    2    3    4    5 
##  6.9 30.6 36.1 16.7  9.7

5.3 Pemodelan

model_ord <- clm(rating ~ temp + contact, data = wine)
summary(model_ord)

## formula: rating ~ temp + contact
## data:    wine
## 
##  link  threshold nobs logLik AIC    niter max.grad cond.H 
##  logit flexible  72   -86.49 184.98 6(0)  4.02e-12 2.7e+01
## 
## Coefficients:
##            Estimate Std. Error z value Pr(>|z|)    
## tempwarm     2.5031     0.5287   4.735 2.19e-06 ***
## contactyes   1.5278     0.4766   3.205  0.00135 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##     Estimate Std. Error z value
## 1|2  -1.3444     0.5171  -2.600
## 2|3   1.2508     0.4379   2.857
## 3|4   3.4669     0.5978   5.800
## 4|5   5.0064     0.7309   6.850

5.4 Odds Ratio dan Interpretasi

or_ord <- data.frame(
  Variabel   = names(coef(model_ord)),
  Koefisien  = round(coef(model_ord), 4),
  Odds_Ratio = round(exp(coef(model_ord)), 4)
)
kable(or_ord, caption = "Koefisien dan Cumulative Odds Ratio — Model Ordinal")

Koefisien dan Cumulative Odds Ratio — Model Ordinal
	Variabel	Koefisien	Odds_Ratio
1\|2	1\|2	-1.3444	0.2607
2\|3	2\|3	1.2508	3.4932
3\|4	3\|4	3.4669	32.0369
4\|5	4\|5	5.0064	149.3667
tempwarm	tempwarm	2.5031	12.2203
contactyes	contactyes	1.5278	4.6080

Interpretasi:

tempwarm: Suhu warm dibandingkan suhu cold — jika koefisiennya positif dan OR > 1, artinya suhu lebih hangat meningkatkan odds anggur dinilai dengan rating kepahitan lebih tinggi.
contactyes: Kontak kulit anggur — jika koefisiennya positif dan OR > 1, artinya kontak kulit meningkatkan odds rating kepahitan yang lebih tinggi.
Threshold (\(\alpha_j\)): Mewakili log-odds kumulatif saat semua prediktor bernilai nol. Nilai threshold meningkat seiring bertambahnya indeks \(j\).

5.5 Uji Asumsi Proportional Odds

model_ord_nom <- clm(rating ~ temp + contact,
                     data = wine, nominal = ~ temp)
anova(model_ord, model_ord_nom)

## Likelihood ratio tests of cumulative link models:
##  
##               formula:                nominal: link: threshold:
## model_ord     rating ~ temp + contact ~1       logit flexible  
## model_ord_nom rating ~ temp + contact ~temp    logit flexible  
## 
##               no.par    AIC  logLik LR.stat df Pr(>Chisq)
## model_ord          6 184.98 -86.492                      
## model_ord_nom      9 187.81 -84.904   3.175  3     0.3654

Interpretasi: Jika \(p\)-value dari likelihood ratio test ini tidak signifikan (\(p > 0.05\)), asumsi proportional odds untuk variabel temp tidak dilanggar — model ordinal yang lebih sederhana sudah memadai. Jika signifikan, diperlukan model yang lebih fleksibel (misalnya model nominal parsial).

5.6 Visualisasi Prediksi per Kombinasi Prediktor

new_data_ord <- expand.grid(
  temp    = c("cold", "warm"),
  contact = c("no", "yes")
)
pred_ord <- predict(model_ord, newdata = new_data_ord, type = "prob")
pred_df  <- cbind(new_data_ord, pred_ord$fit)
colnames(pred_df)[3:7] <- paste0("Rating_", 1:5)

pred_long <- pred_df %>%
  pivot_longer(cols      = starts_with("Rating"),
               names_to  = "Rating",
               values_to = "Probabilitas")

ggplot(pred_long, aes(x = Rating, y = Probabilitas,
                       fill = interaction(temp, contact))) +
  geom_col(position = "dodge", width = 0.7) +
  scale_fill_brewer(palette = "Set2",
                    labels  = c("Cold-No", "Cold-Yes",
                                "Warm-No", "Warm-Yes")) +
  labs(
    title    = "Probabilitas Prediksi Regresi Ordinal per Rating Kepahitan",
    subtitle = "Per kombinasi suhu (cold/warm) dan kontak kulit (no/yes)",
    x    = "Rating Kepahitan",
    y    = "Probabilitas",
    fill = "Temp x Contact"
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

Membaca grafik: Pergeseran distribusi ke kanan (rating lebih tinggi) untuk kondisi warm dan contact = yes mengonfirmasi bahwa kedua prediktor meningkatkan peluang penilaian kepahitan yang lebih tinggi.

6 Regresi Poisson

6.1 Konsep dan Landasan Teori

Regresi Poisson digunakan untuk memodelkan data cacahan (count data) — variabel respons berupa bilangan bulat non-negatif (\(0, 1, 2, \ldots\)), seperti jumlah kejadian, frekuensi kunjungan, atau banyaknya insiden dalam periode tertentu.

Model menggunakan link logaritmik:

\[\log(\lambda_i) = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p \tag{7}\]

sehingga nilai harapan cacahan diprediksi sebagai:

\[\hat\lambda_i = e^{\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p} \tag{8}\]

6.1.1 Rate Ratio

Interpretasi koefisien pada regresi Poisson menggunakan Rate Ratio:

\[RR_k = e^{\beta_k} \tag{9}\]

\(RR = 1\) : tidak ada pengaruh
\(RR > 1\) : meningkatkan laju kejadian
\(RR < 1\) : menurunkan laju kejadian

6.1.2 Asumsi Poisson dan Overdispersi

Asumsi kunci Poisson: \(E(Y) = \text{Var}(Y) = \lambda\) — nilai harapan sama dengan variansi (equidispersion). Jika variansi jauh lebih besar dari nilai harapan, terjadi overdispersi. Dalam kasus ini, model Negative Binomial lebih tepat digunakan.

6.1.3 Notasi

Simbol	Keterangan
\(\lambda_i\)	Laju (rate) kejadian untuk observasi ke-\(i\)
\(e^{\beta_k}\)	Rate Ratio untuk prediktor \(x_k\)
\(\phi\)	Parameter dispersi (Poisson ideal: \(\phi = 1\))

6.2 Dataset: `Affairs` — Jumlah Perselingkuhan

Dataset berisi 601 responden dari survei perilaku pernikahan tahun 1969. Variabel respons affairs adalah jumlah perselingkuhan dalam satu tahun terakhir. Prediktor: age, yearsmarried, religiousness, dan rating (kepuasan pernikahan skala 1–5).

data(Affairs, package = "AER")
kable(head(Affairs, 8), caption = "8 Baris Pertama Dataset Affairs")

8 Baris Pertama Dataset Affairs
	gender	age	yearsmarried	children	religiousness	education	occupation	rating
4	male	37	10.00	no	3	18	7	4
5	female	27	4.00	no	4	14	6	4
11	female	32	15.00	yes	1	12	1	4
16	male	57	15.00	yes	5	18	6	5
23	male	22	0.75	no	2	17	6	3
29	female	32	1.50	no	2	17	5	5
44	female	22	0.75	no	2	12	1	3
45	male	57	15.00	yes	2	14	4	4

cat("Distribusi variabel respons (affairs):\n")

## Distribusi variabel respons (affairs):

print(table(Affairs$affairs))

## 
##   0   1   2   3   7  12 
## 451  34  17  19  42  38

ggplot(Affairs, aes(x = affairs)) +
  geom_histogram(bins = 12, fill = "#2980b9", color = "white", alpha = 0.85) +
  labs(
    title    = "Distribusi Jumlah Perselingkuhan",
    subtitle = "Data sangat condong ke kanan (right-skewed) — khas distribusi Poisson",
    x = "Jumlah Kejadian (affairs)",
    y = "Frekuensi"
  ) +
  theme_minimal(base_size = 13)

6.3 Pemodelan

model_pois <- glm(affairs ~ age + yearsmarried + religiousness + rating,
                  family = poisson(link = "log"),
                  data   = Affairs)
summary(model_pois)

## 
## Call:
## glm(formula = affairs ~ age + yearsmarried + religiousness + 
##     rating, family = poisson(link = "log"), data = Affairs)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    2.748394   0.188189  14.604  < 2e-16 ***
## age           -0.027057   0.005686  -4.759 1.95e-06 ***
## yearsmarried   0.110078   0.009812  11.219  < 2e-16 ***
## religiousness -0.360786   0.030869 -11.688  < 2e-16 ***
## rating        -0.401699   0.027285 -14.722  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2925.5  on 600  degrees of freedom
## Residual deviance: 2377.5  on 596  degrees of freedom
## AIC: 2881.5
## 
## Number of Fisher Scoring iterations: 7

6.4 Rate Ratio dan Confidence Interval

rr_table <- data.frame(
  Variabel   = names(coef(model_pois)),
  Koefisien  = round(coef(model_pois), 4),
  Rate_Ratio = round(exp(coef(model_pois)), 4),
  CI_Lower   = round(exp(confint(model_pois))[, 1], 4),
  CI_Upper   = round(exp(confint(model_pois))[, 2], 4)
)
kable(rr_table, caption = "Rate Ratio (exp(beta)) dan Confidence Interval 95% — Regresi Poisson")

Rate Ratio (exp(beta)) dan Confidence Interval 95% — Regresi Poisson
	Variabel	Koefisien	Rate_Ratio	CI_Lower	CI_Upper
(Intercept)	(Intercept)	2.7484	15.6175	10.7978	22.5820
age	age	-0.0271	0.9733	0.9624	0.9841
yearsmarried	yearsmarried	0.1101	1.1164	1.0951	1.1381
religiousness	religiousness	-0.3608	0.6971	0.6561	0.7405
rating	rating	-0.4017	0.6692	0.6344	0.7060

Interpretasi:

religiousness: Setiap kenaikan 1 unit tingkat religiusitas, laju perselingkuhan berubah sebesar faktor \(e^{\hat\beta}\). Jika RR < 1 dan signifikan, religiusitas yang lebih tinggi menurunkan laju perselingkuhan.
rating: Kepuasan pernikahan yang lebih tinggi diharapkan memiliki RR < 1 — semakin puas dengan pernikahan, semakin rendah laju perselingkuhan.
yearsmarried: Lama pernikahan yang lebih panjang dapat meningkatkan laju, sejalan dengan teori bahwa kejenuhan bertambah seiring waktu.
age: Pengaruh usia bisa berkorelasi dengan yearsmarried — perlu hati-hati dalam interpretasi akibat potensi multikolinieritas.

6.5 Uji Overdispersi

dispersiontest(model_pois)

## 
##  Overdispersion test
## 
## data:  model_pois
## z = 5.4005, p-value = 3.322e-08
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   6.816143

model_nb <- glm.nb(affairs ~ age + yearsmarried + religiousness + rating,
                   data = Affairs)

cat("=== Perbandingan AIC ===\n")

## === Perbandingan AIC ===

cat("Poisson AIC        :", round(AIC(model_pois), 2), "\n")

## Poisson AIC        : 2881.53

cat("Neg. Binomial AIC  :", round(AIC(model_nb),   2), "\n")

## Neg. Binomial AIC  : 1469.22

if (AIC(model_nb) < AIC(model_pois)) {
  cat("\n>> Model Negative Binomial lebih baik (AIC lebih kecil).\n")
  cat(">> Ini mengindikasikan adanya OVERDISPERSI pada data.\n")
} else {
  cat("\n>> Model Poisson sudah cukup memadai.\n")
}

## 
## >> Model Negative Binomial lebih baik (AIC lebih kecil).
## >> Ini mengindikasikan adanya OVERDISPERSI pada data.

Tindak lanjut: Jika dispersiontest signifikan (\(p < 0.05\)) dan AIC Negative Binomial lebih kecil, gunakan model Negative Binomial sebagai model final karena asumsi equidispersi Poisson dilanggar.

6.6 Visualisasi: Efek Rating Pernikahan

eff_df <- as.data.frame(effect("rating", model_pois,
                                xlevels = list(rating = 1:5)))

ggplot(eff_df, aes(x = rating, y = fit)) +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              fill = "#3498db", alpha = 0.25) +
  geom_line(color  = "#2980b9", linewidth = 1.4) +
  geom_point(color = "#1a5276", size = 3.5) +
  labs(
    title    = "Efek Rating Pernikahan terhadap Laju Perselingkuhan",
    subtitle = "Regresi Poisson — variabel lain di-hold pada nilai rata-rata",
    x = "Rating Pernikahan (1 = sangat tidak bahagia, 5 = sangat bahagia)",
    y = "Prediksi Laju Perselingkuhan (lambda)"
  ) +
  theme_minimal(base_size = 13)

Membaca grafik: Kurva yang menurun dari rating 1 ke 5 mengonfirmasi hubungan negatif antara kepuasan pernikahan dan laju perselingkuhan. Pita biru mewakili interval kepercayaan 95%.

7 Ringkasan dan Perbandingan Model

ringkasan <- data.frame(
  Jenis_Regresi = c("Biner", "Multinomial", "Ordinal", "Poisson"),
  Respons       = c("Dikotom (0/1)", "Nominal >= 3 kat.", "Ordinal >= 3 kat.", "Count (cacahan)"),
  Link_Function = c("Logit", "Softmax/Logit", "Cumulative Logit", "Log"),
  Interpretasi  = c("Odds Ratio", "OR vs referensi", "Cum. Odds Ratio", "Rate Ratio"),
  Package       = c("stats (base R)", "nnet", "ordinal", "stats + AER"),
  Dataset       = c("Pima.tr", "hsbdemo (UCLA)", "wine", "Affairs"),
  AIC           = c(
    round(AIC(model_bin),  2),
    round(AIC(model_mn),   2),
    round(AIC(model_ord),  2),
    round(AIC(model_pois), 2)
  )
)
kable(ringkasan, caption = "Ringkasan Keempat Model Regresi yang Dianalisis")

Ringkasan Keempat Model Regresi yang Dianalisis
Jenis_Regresi	Respons	Link_Function	Interpretasi	Package	Dataset	AIC
Biner	Dikotom (0/1)	Logit	Odds Ratio	stats (base R)	Pima.tr	198.27
Multinomial	Nominal >= 3 kat.	Softmax/Logit	OR vs referensi	nnet	hsbdemo (UCLA)	370.34
Ordinal	Ordinal >= 3 kat.	Cumulative Logit	Cum. Odds Ratio	ordinal	wine	184.98
Poisson	Count (cacahan)	Log	Rate Ratio	stats + AER	Affairs	2881.53

Catatan AIC: AIC hanya bermakna untuk membandingkan model pada dataset yang sama. Nilai AIC lintas dataset berbeda tidak bisa dibandingkan secara langsung — masing-masing hanya mencerminkan goodness-of-fit relatif dalam konteksnya sendiri.

8 Panduan Pemilihan Model

Kondisi Variabel Respons	Model yang Tepat
2 kategori (ya/tidak, sukses/gagal)	Regresi Logistik Biner
>= 3 kategori tanpa urutan (merah/biru/hijau)	Regresi Multinomial
>= 3 kategori dengan urutan (rendah/sedang/tinggi)	Regresi Ordinal
Bilangan bulat >= 0 (0, 1, 2, 3, …)	Regresi Poisson
Count dengan variansi >> mean	Negative Binomial (alternatif Poisson)

9 Referensi

Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Wiley.
Venables, W. N. & Ripley, B. D. (2002). Modern Applied Statistics with S. Springer. — Package MASS
Ripley, B. (2023). nnet: Feed-Forward Neural Networks and Multinomial Log-Linear Models. — Package nnet
Christensen, R. H. B. (2023). ordinal: Regression Models for Ordinal Data. — Package ordinal
Kleiber, C. & Zeileis, A. (2008). Applied Econometrics with R. Springer. — Package AER
UCLA OARC Statistics. R Data Analysis Examples. https://stats.oarc.ucla.edu/r/

Laporan Praktikum Pemodelan Regresi

Regresi Logistik Biner | Multinomial | Ordinal | Poisson

Bijan Ramadhan Aditri

01 June 2026

1 Pendahuluan

2 Instalasi dan Pemuatan Package

3 Regresi Logistik Biner

3.1 Konsep dan Landasan Teori

3.1.1 Notasi

3.1.2 Interpretasi Odds Ratio

3.2 Dataset: Pima.tr — Diabetes Suku Pima

3.3 Pemodelan

3.4 Odds Ratio dan Confidence Interval

3.5 Visualisasi: Kurva Logistik

3.6 Evaluasi Model

4 Regresi Multinomial

4.1 Konsep dan Landasan Teori

4.1.1 Notasi

4.2 Dataset: hsbdemo — Program Pendidikan Siswa

4.3 Pemodelan

4.4 Koefisien, Z-Score, dan Odds Ratio

4.5 Visualisasi Prediksi

5 Regresi Ordinal

5.1 Konsep dan Landasan Teori

5.1.1 Asumsi Proportional Odds

5.1.2 Notasi

5.2 Dataset: wine — Penilaian Kepahitan Anggur

5.3 Pemodelan

5.4 Odds Ratio dan Interpretasi

5.5 Uji Asumsi Proportional Odds

5.6 Visualisasi Prediksi per Kombinasi Prediktor

6 Regresi Poisson

6.1 Konsep dan Landasan Teori

6.1.1 Rate Ratio

6.1.2 Asumsi Poisson dan Overdispersi

6.1.3 Notasi

6.2 Dataset: Affairs — Jumlah Perselingkuhan

6.3 Pemodelan

6.4 Rate Ratio dan Confidence Interval

6.5 Uji Overdispersi

6.6 Visualisasi: Efek Rating Pernikahan

7 Ringkasan dan Perbandingan Model

8 Panduan Pemilihan Model

9 Referensi

3.2 Dataset: `Pima.tr` — Diabetes Suku Pima

4.2 Dataset: `hsbdemo` — Program Pendidikan Siswa

5.2 Dataset: `wine` — Penilaian Kepahitan Anggur

6.2 Dataset: `Affairs` — Jumlah Perselingkuhan