Analisis Regresi Logistik: Biner, Multinomial, Ordinal, dan Poisson

---

1 Regresi Logistik Biner

Regresi logistik biner adalah model regresi untuk respons kategorik dengan tepat dua kategori, yaitu kejadian yang diamati terjadi (\(Y = 1\)) atau tidak terjadi (\(Y = 0\)).

Misalkan respons: \[Y_i \in \{0, 1\}\] Model memodelkan log-odds (logit) dari peluang kejadian sebagai fungsi linear prediktor:

\[ \log\left(\frac{P(Y_i = 1 \mid \mathbf{x}_i)}{P(Y_i = 0 \mid \mathbf{x}_i)}\right) = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}. \]

Dari persamaan ini, peluang kejadian diperoleh melalui fungsi logistik:

\[ P(Y_i = 1 \mid \mathbf{x}_i) = \frac{\exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip})} {1 + \exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip})}. \]

Artinya, model membangun satu persamaan logit yang membandingkan peluang kejadian (\(Y = 1\)) terhadap peluang tidak terjadi (\(Y = 0\)) sebagai kategori referensi.

Studi kasus: Analisis sikap publik terhadap hukuman mati berdasarkan karakteristik demografis dan pandangan politik responden.

Tujuan analisis: Membangun model klasifikasi untuk memprediksi apakah seorang responden mendukung (FAVOR) atau menolak (OPPOSE) hukuman mati berdasarkan usia, jenis kelamin, dan pandangan politik.

1.1 Ringkasan Masalah

Hukuman mati merupakan isu sosial yang sangat kontroversial dan banyak diteliti dalam ilmu sosial. Preferensi seseorang terhadap hukuman mati sering kali dikaitkan dengan latar belakang demografis dan orientasi politik. Kasus ini cocok dianalisis dengan regresi logistik biner karena variabel respons bersifat biner:

• 1 = FAVOR (mendukung hukuman mati)
• 0 = OPPOSE (menolak hukuman mati)

Model yang digunakan adalah:

\[ \log\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots + \beta_k X_{ki}, \]

dengan \(p_i = P(Y_i=1 \mid X_i)\) adalah peluang responden ke-\(i\) mendukung hukuman mati (FAVOR).

1.2 Sumber Data

Data yang digunakan dalam analisis ini bersumber dari General Social Survey (GSS) tahun 2024, sebuah survei opini publik berskala nasional di Amerika Serikat (United States) yang mengukur tren sosial, sikap, dan karakteristik demografis masyarakat. Fokus amatan pada penelitian ini mencakup sikap responden terhadap hukuman mati beserta karakteristik demografis dan pandangan politik mereka. Dataset yang dianalisis terdiri dari 591 responden dengan 5 variabel.

1.3 Persiapan Data

1.3.1 Membaca Data dari Excel

raw_data <- readxl::read_excel("Data Biner.xlsx", sheet = 1)

ringkasan_data <- data.frame(
  Keterangan = c("Jumlah observasi", "Jumlah variabel"),
  Nilai      = c(nrow(raw_data), ncol(raw_data))
)

knitr::kable(
  ringkasan_data,
  caption = "Ukuran dataset"
)

Ukuran dataset

Keterangan	Nilai
Jumlah observasi	591
Jumlah variabel	5

Interpretasi output: Dataset berisi 591 observasi dan 5 variabel. Jumlah observasi ini memadai untuk membangun model regresi logistik biner dengan tiga prediktor. Variabel respons adalah cappun, yaitu sikap responden terhadap hukuman mati (FAVOR atau OPPOSE).

1.3.2 Penjelasan Variabel

kamus_variabel <- data.frame(
  `Nama kolom`          = c("id_", "age", "sex", "polviews", "cappun"),
  `Keterangan`          = c(
    "Nomor identifikasi responden",
    "Usia responden (tahun)",
    "Jenis kelamin responden",
    "Pandangan politik responden",
    "Sikap terhadap hukuman mati (variabel respons)"
  ),
  `Kategori / Rentang`  = c(
    "—",
    "18 – 89 tahun",
    "MALE, FEMALE",
    "Conservative, Liberal",
    "FAVOR (mendukung), OPPOSE (menolak)"
  ),
  `Tipe dalam analisis` = c(
    "ID (tidak dipakai dalam model)",
    "Numerik",
    "Kategorik",
    "Kategorik",
    "Respons biner"
  ),
  check.names = FALSE
)

knitr::kable(
  kamus_variabel,
  caption = "Penjelasan variabel dataset sikap terhadap hukuman mati"
)

Penjelasan variabel dataset sikap terhadap hukuman mati

Nama kolom	Keterangan	Kategori / Rentang	Tipe dalam analisis
id_	Nomor identifikasi responden	—	ID (tidak dipakai dalam model)
age	Usia responden (tahun)	18 – 89 tahun	Numerik
sex	Jenis kelamin responden	MALE, FEMALE	Kategorik
polviews	Pandangan politik responden	Conservative, Liberal	Kategorik
cappun	Sikap terhadap hukuman mati (variabel respons)	FAVOR (mendukung), OPPOSE (menolak)	Respons biner

Interpretasi output: Tabel penjelasan variabel menjelaskan makna setiap kolom. Variabel age merupakan satu-satunya prediktor numerik. Variabel sex dan polviews adalah prediktor kategorik, sedangkan cappun adalah variabel respons biner yang menjadi target prediksi.

1.3.3 Transformasi dan Pembersihan Data

library(dplyr)
data_bersih <- raw_data |>
  mutate(
    # Usia sebagai numerik
    age      = as.numeric(age), 
    # Jenis kelamin sebagai factor; MALE sebagai referensi
    sex = factor(sex, levels = c("MALE", "FEMALE")),
    # Pandangan politik sebagai factor; Conservative sebagai referensi
    polviews = factor(polviews, levels = c("Conservative", "Liberal")),
    # Respons biner: 1 = FAVOR, 0 = OPPOSE
    Y_biner = as.integer(cappun == "FAVOR"),
    # Label factor untuk visualisasi
    Y_label = factor(
      ifelse(Y_biner == 1, "FAVOR (Mendukung)", "OPPOSE (Menolak)"),
      levels = c("OPPOSE (Menolak)", "FAVOR (Mendukung)")
    )
  ) |>
  tidyr::drop_na()

knitr::kable(
  head(data_bersih %>%
    select(id_, age, sex, polviews, cappun, Y_biner), 8),
  caption = "Contoh data setelah transformasi"
)

Contoh data setelah transformasi

id_	age	sex	polviews	cappun	Y_biner
7	48	FEMALE	Liberal	OPPOSE	0
11	23	FEMALE	Liberal	FAVOR	1
24	48	MALE	Conservative	FAVOR	1
37	60	FEMALE	Conservative	FAVOR	1
45	58	MALE	Conservative	FAVOR	1
50	65	MALE	Conservative	FAVOR	1
60	63	FEMALE	Liberal	FAVOR	1
61	57	FEMALE	Liberal	OPPOSE	0

Interpretasi output: Tabel menampilkan delapan baris pertama data setelah transformasi. Kolom Y_biner menunjukkan nilai 1 untuk responden yang mendukung (FAVOR) dan 0 untuk yang menolak (OPPOSE) hukuman mati. Kategori referensi yang dipilih adalah MALE untuk sex dan Conservative untuk polviews, sehingga odds ratio untuk level lain dibandingkan terhadap kategori-kategori ini.

1.3.4 Distribusi Kelas Respons

class_summary <- data_bersih %>%
  count(Y_label, name = "Jumlah") %>%
  mutate(
    Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1)
  ) %>%
  rename(Status = Y_label)

knitr::kable(
  class_summary,
  caption = "Distribusi kelas respons"
)

Distribusi kelas respons

Status	Jumlah	Proporsi
OPPOSE (Menolak)	231	39.1%
FAVOR (Mendukung)	360	60.9%

Interpretasi output: Dari 591 responden, sebanyak 360 (60,9%) mendukung hukuman mati dan 231 (39,1%) menolaknya. Kelas respons tidak seimbang sempurna, tetapi proporsi kelas minoritas (OPPOSE) masih cukup besar sehingga model klasifikasi dapat mempelajari pola kedua kelompok dengan baik.

1.4 Eksplorasi Singkat

1.4.1 Distribusi Variabel Respons

ggplot(data_bersih, aes(x = Y_label, fill = Y_label)) +
  geom_bar(width = 0.62, color = "white", linewidth = 0.8) +
  geom_text(
    stat     = "count",
    aes(label = after_stat(count)),
    vjust    = -0.4,
    fontface = "bold"
  ) +
  scale_fill_manual(
    values = c("OPPOSE (Menolak)" = "#2a9d8f",
               "FAVOR (Mendukung)" = "#e76f51")
  ) +
  labs(
    title = "Distribusi Sikap terhadap Hukuman Mati",
    x     = NULL,
    y     = "Jumlah responden"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

Interpretasi output: Grafik batang memperlihatkan bahwa responden yang mendukung hukuman mati (FAVOR) lebih banyak daripada yang menolak (OPPOSE). Ini mencerminkan pola umum dalam survei opini publik di mana dukungan terhadap hukuman mati seringkali merupakan posisi mayoritas, terutama pada kelompok yang lebih konservatif dan lebih tua.

1.4.2 Ringkasan Variabel Numerik (Usia)

numeric_summary <- data_bersih %>%
  group_by(Y_label) %>%
  summarise(
    Jumlah           = n(),
    Ratarata_usia = mean(age, na.rm = TRUE),
    Median_usia   = median(age, na.rm = TRUE),
    Std_deviasi   = sd(age, na.rm = TRUE),
    Min_usia      = min(age),
    Max_usia       = max(age),
    .groups = "drop"
  ) %>%
  rename(Status = Y_label) %>%
  mutate(across(where(is.numeric), round, 2))

knitr::kable(
  numeric_summary,
  caption = "Ringkasan usia responden menurut sikap terhadap hukuman mati"
)

Ringkasan usia responden menurut sikap terhadap hukuman mati

Status	Jumlah	Ratarata_usia	Median_usia	Std_deviasi	Min_usia	Max_usia
OPPOSE (Menolak)	231	51.68	55	18.98	18	89
FAVOR (Mendukung)	360	53.44	55	17.00	18	89

Interpretasi output: Responden yang mendukung hukuman mati (FAVOR) memiliki rata-rata usia lebih tinggi dibandingkan yang menolak (OPPOSE). Perbedaan rata-rata usia ini memberikan indikasi awal bahwa usia berpotensi menjadi prediktor yang relevan dalam model, di mana responden yang lebih tua cenderung lebih mendukung hukuman mati.

1.4.3 Distribusi Usia dan Pandangan Politik

ggplot(
  data_bersih,
  aes(x = age, fill = Y_label)
) +
  geom_density(alpha = 0.55, color = "white", linewidth = 0.7) +
  facet_wrap(~ polviews, labeller = label_both) +
  scale_fill_manual(
    values = c("OPPOSE (Menolak)" = "#2a9d8f",
               "FAVOR (Mendukung)" = "#e76f51")
  ) +
  labs(
    title    = "Distribusi Usia Responden Menurut Pandangan Politik dan Sikap",
    subtitle = "Dikelompokkan berdasarkan pandangan politik",
    x        = "Usia (tahun)",
    y        = "Kepadatan",
    fill     = "Sikap"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Interpretasi output: Grafik kepadatan yang difacet berdasarkan pandangan politik menunjukkan pola distribusi usia untuk kelompok FAVOR dan OPPOSE pada responden Conservative maupun Liberal. Jika dua kurva terpisah jauh dalam suatu facet, artinya usia menjadi pembeda yang kuat dalam kelompok tersebut. Perbedaan pola antara facet Conservative dan Liberal juga menggambarkan peran interaksi antara usia dan pandangan politik.

1.4.4 Proporsi Dukungan Berdasarkan Jenis Kelamin dan Pandangan Politik

prop_summary <- data_bersih %>%
  group_by(sex, polviews) %>%
  summarise(
    prop_favor = mean(Y_biner),
    n          = n(),
    .groups    = "drop"
  ) %>%
  mutate(
    label = paste0(round(prop_favor * 100, 1), "%\n(n=", n, ")")
  )

ggplot(prop_summary,
       aes(x = sex, y = prop_favor, fill = polviews)) +
  geom_col(position = "dodge", width = 0.65,
           color = "white", linewidth = 0.8) +
  geom_text(
    aes(label = label),
    position = position_dodge(width = 0.65),
    vjust    = -0.3,
    size     = 3.2
  ) +
  scale_fill_manual(values = c("Conservative" = "#e76f51",
                               "Liberal"       = "#2a9d8f")) +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1.1)) +
  labs(
    title    = "Proporsi Dukungan Hukuman Mati (FAVOR)",
    subtitle = "Menurut jenis kelamin dan pandangan politik",
    x        = "Jenis kelamin",
    y        = "Proporsi mendukung (FAVOR)",
    fill     = "Pandangan politik"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Interpretasi output: Grafik batang berkelompok memperlihatkan proporsi dukungan hukuman mati untuk setiap kombinasi jenis kelamin dan pandangan politik. Responden Conservative umumnya memiliki proporsi FAVOR lebih tinggi dibandingkan Liberal pada kelompok jenis kelamin yang sama. Perbedaan proporsi ini mengindikasikan bahwa polviews merupakan prediktor yang penting secara substantif.

1.5 Pembagian Training dan Testing

Data dibagi menjadi 80% training dan 20% testing secara stratified agar proporsi FAVOR dan OPPOSE tetap seimbang pada kedua subset.

stratified_split <- function(y, prop = 0.8) {
  idx_by_class <- split(seq_along(y), y)
  train_idx <- lapply(
    idx_by_class,
    function(idx) sample(idx, size = floor(length(idx) * prop))
  )
  unlist(train_idx, use.names = FALSE)
}

set.seed(42)
train_id   <- stratified_split(data_bersih$Y_biner, prop = 0.8)
train_data <- data_bersih[ train_id, ]
test_data  <- data_bersih[-train_id, ]

split_summary <- bind_rows(
  train_data %>% count(Y_label) %>% mutate(data = "Training"),
  test_data  %>% count(Y_label) %>% mutate(data = "Testing")
) %>%
  group_by(data) %>%
  mutate(proporsi = scales::percent(n / sum(n), accuracy = 0.1)) %>%
  ungroup() %>%
  select(data, Y_label, n, proporsi) %>%
  rename(
    Data     = data,
    Status   = Y_label,
    Jumlah   = n,
    Proporsi = proporsi
  )

knitr::kable(split_summary,
             caption = "Distribusi kelas pada training dan testing")

Distribusi kelas pada training dan testing

Data	Status	Jumlah	Proporsi
Training	OPPOSE (Menolak)	184	39.0%
Training	FAVOR (Mendukung)	288	61.0%
Testing	OPPOSE (Menolak)	47	39.5%
Testing	FAVOR (Mendukung)	72	60.5%

Interpretasi output: Data training berisi 472 observasi dan data testing berisi 119 observasi. Karena split dilakukan secara stratified, proporsi FAVOR dan OPPOSE pada training maupun testing tetap mendekati proporsi data awal, sehingga evaluasi pada testing tidak bias akibat perbedaan komposisi kelas.

1.6 Model Regresi Logistik Biner

1.6.1 Rumus Model Logistik

Misalkan \(Y_i\) adalah sikap responden ke-\(i\) terhadap hukuman mati, dengan:

\[ Y_i = \begin{cases} 1, & \text{jika mendukung hukuman mati (FAVOR)},\\ 0, & \text{jika menolak hukuman mati (OPPOSE)}. \end{cases} \]

Peluang mendukung hukuman mati dinotasikan sebagai:

\[ p_i = P(Y_i = 1 \mid X_i). \]

Pada regresi logistik biner, peluang tersebut dimodelkan melalui fungsi logit:

\[ \text{logit}(p_i) = \ln\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1 \cdot \text{age}_i + \beta_2 \cdot \mathbf{1}[\text{sex}_i = \text{FEMALE}] + \beta_3 \cdot \mathbf{1}[\text{polviews}_i = \text{Liberal}]. \]

Bentuk peluang prediksi diperoleh dengan transformasi balik:

\[ \hat{p}_i = \frac{\exp(\hat{\eta}_i)}{1+\exp(\hat{\eta}_i)} = \frac{1}{1+\exp(-\hat{\eta}_i)}. \]

Interpretasi odds ratio: Untuk prediktor numerik \(X_j\), \(OR_j = \exp(\beta_j)\). Jika \(OR_j > 1\), kenaikan prediktor menaikkan odds mendukung hukuman mati; jika \(OR_j < 1\), sebaliknya. Untuk prediktor kategorik, odds ratio dibandingkan terhadap kategori referensi (MALE untuk sex, Conservative untuk polviews).

1.6.2 Estimasi Model pada Data Training

model_fit <- glm(
  Y_biner ~ age + sex + polviews,
  data   = train_data,
  family = binomial(link = "logit")
)

ringkasan_model <- data.frame(
  Keterangan = c(
    "Jumlah observasi training",
    "Null deviance",
    "Residual deviance",
    "Derajat bebas residual",
    "AIC"
  ),
  Nilai = c(
    nobs(model_fit),
    round(model_fit$null.deviance, 3),
    round(model_fit$deviance, 3),
    model_fit$df.residual,
    round(AIC(model_fit), 3)
  )
)

knitr::kable(
  ringkasan_model,
  caption = "Ringkasan kecocokan model regresi logistik"
)

Ringkasan kecocokan model regresi logistik

Keterangan	Nilai
Jumlah observasi training	472.000
Null deviance	631.227
Residual deviance	563.852
Derajat bebas residual	468.000
AIC	571.852

Interpretasi output: Model diestimasi menggunakan 472 observasi training. Nilai residual deviance (563.852) lebih kecil daripada null deviance (631.227), yang berarti ketiga prediktor (usia, jenis kelamin, pandangan politik) secara bersama-sama memberikan perbaikan kecocokan dibandingkan model kosong yang hanya memakai intercept. Nilai AIC sebesar 571.852 dapat digunakan untuk membandingkan model ini dengan model alternatif.

1.6.3 Odds Ratio

coef_table <- broom::tidy(model_fit) %>%
  mutate(
    odds_ratio = exp(estimate),
    ci_low     = exp(estimate - 1.96 * std.error),
    ci_high    = exp(estimate + 1.96 * std.error)
  ) %>%
  arrange(p.value) %>%
  transmute(
    `Variabel/level`                 = term,
    `Odds ratio`                     = round(odds_ratio, 3),
    `Interval kepercayaan 95 persen` = paste0(round(ci_low, 3),
                                              " - ",
                                              round(ci_high, 3)),
    `p-value`                        = signif(p.value, 3)
  )

knitr::kable(
  coef_table,
  caption = "Koefisien model: odds ratio dan interval kepercayaan 95%"
)

Koefisien model: odds ratio dan interval kepercayaan 95%

Variabel/level	Odds ratio	Interval kepercayaan 95 persen	p-value
polviewsLiberal	0.202	0.134 - 0.303	0.00e+00
(Intercept)	5.369	2.523 - 11.425	1.29e-05
sexFEMALE	0.703	0.471 - 1.05	8.53e-02
age	0.995	0.984 - 1.007	4.26e-01

Interpretasi output: Tabel odds ratio diurutkan dari koefisien dengan p-value terkecil.

• polviewsLiberal: OR < 1 dan p-value kecil menunjukkan bahwa responden dengan pandangan Liberal memiliki odds mendukung hukuman mati yang secara signifikan lebih rendah dibandingkan responden Conservative (kategori referensi). Ini merupakan prediktor terkuat dalam model.
• sexFEMALE: jika OR < 1, responden perempuan memiliki odds mendukung hukuman mati lebih rendah dibandingkan laki-laki, atau sebaliknya jika OR > 1.
• age: OR mendekati 1 mengindikasikan efek usia per tahun. Meski kecil, efek kumulatif selama beberapa dekade bisa signifikan secara substantif.
• Interval kepercayaan yang tidak melewati nilai 1 dan p-value < 0,05 memberikan bukti statistik yang kuat bahwa prediktor tersebut relevan dalam model.

1.7 Prediksi dan Evaluasi Klasifikasi

1.7.1 Rumus Prediksi Kelas

Setelah model menghasilkan peluang prediksi \(\hat{p}_i\), kelas ditentukan menggunakan threshold \(c\):

\[ \hat{Y}_i = \begin{cases} 1, & \text{jika } \hat{p}_i \ge c,\\ 0, & \text{jika } \hat{p}_i < c. \end{cases} \]

1.7.2 Rumus Confusion Matrix dan Metrik Evaluasi

Notasi confusion matrix yang digunakan:

• TP / true positive: aktual FAVOR dan diprediksi FAVOR.
• TN / true negative: aktual OPPOSE dan diprediksi OPPOSE.
• FP / false positive: aktual OPPOSE tetapi diprediksi FAVOR.
• FN / false negative: aktual FAVOR tetapi diprediksi OPPOSE.

\[ \text{Akurasi} = \frac{TP + TN}{TP + TN + FP + FN} \qquad \text{Error rate} = 1 - \text{Akurasi} \]

\[ \text{Sensitivity} = \frac{TP}{TP + FN} \qquad \text{Specificity} = \frac{TN}{TN + FP} \qquad \text{Precision} = \frac{TP}{TP + FP} \]

\[ \text{F1-score} = 2 \times \frac{\text{Precision} \times \text{Sensitivity}} {\text{Precision} + \text{Sensitivity}} \qquad \text{Balanced accuracy} = \frac{\text{Sensitivity} + \text{Specificity}}{2} \]

Dalam konteks ini, sensitivity mengukur kemampuan model mendeteksi responden yang benar-benar mendukung hukuman mati (FAVOR), sedangkan specificity mengukur kemampuan model mengenali responden yang menolak (OPPOSE).

safe_div <- function(num, den) {
  ifelse(den == 0, NA_real_, num / den)
}

classification_metrics <- function(actual, prob, threshold = 0.5) {
  pred <- as.integer(prob >= threshold)
  tp   <- sum(pred == 1 & actual == 1)
  tn   <- sum(pred == 0 & actual == 0)
  fp   <- sum(pred == 1 & actual == 0)
  fn   <- sum(pred == 0 & actual == 1)

  sensitivity <- safe_div(tp, tp + fn)
  specificity <- safe_div(tn, tn + fp)
  precision   <- safe_div(tp, tp + fp)
  npv         <- safe_div(tn, tn + fn)
  accuracy    <- safe_div(tp + tn, tp + tn + fp + fn)

  data.frame(
    threshold                 = threshold,
    accuracy                  = accuracy,
    error_rate                = 1 - accuracy,
    sensitivity               = sensitivity,
    specificity               = specificity,
    precision                 = precision,
    negative_predictive_value = npv,
    f1_score                  = safe_div(2 * precision * sensitivity,
                                         precision + sensitivity),
    balanced_accuracy         = (sensitivity + specificity) / 2,
    false_positive_rate       = 1 - specificity,
    false_negative_rate       = 1 - sensitivity
  )
}

format_metrics_indonesia <- function(x) {
  x %>%
    transmute(
      Threshold           = threshold,
      Akurasi             = accuracy,
      `Error rate`        = error_rate,
      Sensitivity         = sensitivity,
      Specificity         = specificity,
      Presisi             = precision,
      NPV                 = negative_predictive_value,
      `F1-score`          = f1_score,
      `Balanced accuracy` = balanced_accuracy,
      FPR                 = false_positive_rate,
      FNR                 = false_negative_rate
    )
}

confusion_matrix <- function(actual, prob, threshold = 0.5) {
  pred_label <- factor(
    ifelse(prob >= threshold, "Prediksi FAVOR", "Prediksi OPPOSE"),
    levels = c("Prediksi FAVOR", "Prediksi OPPOSE")
  )
  actual_label <- factor(
    ifelse(actual == 1, "Aktual FAVOR", "Aktual OPPOSE"),
    levels = c("Aktual FAVOR", "Aktual OPPOSE")
  )
  addmargins(table(actual_label, pred_label))
}

1.7.3 Prediksi Peluang

p_train <- predict(model_fit, newdata = train_data, type = "response")
p_test  <- predict(model_fit, newdata = test_data,  type = "response")

prediction_preview <- head(
  data.frame(
    `Status aktual`           = test_data$Y_label,
    `Usia`                    = test_data$age,
    `Jenis kelamin`           = test_data$sex,
    `Pandangan politik`       = test_data$polviews,
    `Peluang prediksi FAVOR`  = round(p_test, 4),
    check.names = FALSE
  ),
  8
)

knitr::kable(
  prediction_preview,
  caption = "Contoh peluang prediksi mendukung hukuman mati (FAVOR) pada data testing"
)

Contoh peluang prediksi mendukung hukuman mati (FAVOR) pada data testing

Status aktual	Usia	Jenis kelamin	Pandangan politik	Peluang prediksi FAVOR
FAVOR (Mendukung)	58	MALE	Conservative	0.8034
OPPOSE (Menolak)	79	MALE	Liberal	0.4274
FAVOR (Mendukung)	63	FEMALE	Liberal	0.3614
OPPOSE (Menolak)	55	MALE	Liberal	0.4552
FAVOR (Mendukung)	36	MALE	Liberal	0.4775
OPPOSE (Menolak)	72	MALE	Liberal	0.4355
OPPOSE (Menolak)	80	MALE	Liberal	0.4263
OPPOSE (Menolak)	21	MALE	Conservative	0.8295

Interpretasi output: Tabel ini menampilkan beberapa contoh observasi testing beserta peluang prediksi mendukung hukuman mati. Semakin besar nilai peluang, semakin tinggi kecenderungan responden diklasifikasikan sebagai FAVOR. Peluang ini belum menjadi kelas akhir sampai dibandingkan dengan threshold.

1.7.4 Evaluasi pada Threshold 0.50

confusion_default <- confusion_matrix(test_data$Y_biner, p_test,
                                      threshold = 0.5)
metrics_default   <- classification_metrics(test_data$Y_biner, p_test,
                                            threshold = 0.5) %>%
  format_metrics_indonesia() %>%
  mutate(across(where(is.numeric), round, 3))

knitr::kable(
  confusion_default,
  caption = "Confusion matrix testing pada threshold 0.50"
)

Confusion matrix testing pada threshold 0.50

	Prediksi FAVOR	Prediksi OPPOSE	Sum
Aktual FAVOR	49	23	72
Aktual OPPOSE	11	36	47
Sum	60	59	119

knitr::kable(
  metrics_default,
  caption = "Metrik testing pada threshold 0.50"
)

Metrik testing pada threshold 0.50

Threshold	Akurasi	Error rate	Sensitivity	Specificity	Presisi	NPV	F1-score	Balanced accuracy	FPR	FNR
0.5	0.714	0.286	0.681	0.766	0.817	0.61	0.742	0.723	0.234	0.319

Interpretasi output: Pada threshold 0.50, confusion matrix menunjukkan berapa responden FAVOR yang berhasil dideteksi dan berapa responden OPPOSE yang tetap dikenali sebagai OPPOSE. Akurasi menyatakan proporsi total prediksi yang benar. Jika sensitivity masih rendah, model melewatkan banyak responden yang sebenarnya mendukung hukuman mati. Jika specificity rendah, banyak responden yang sebenarnya menolak justru diklasifikasikan sebagai pendukung.

1.8 Kurva ROC dan Threshold Optimal

Setiap titik pada kurva ROC berasal dari satu nilai threshold \(c\):

\[ TPR(c) = \text{Sensitivity}(c) = \frac{TP(c)}{TP(c)+FN(c)} \qquad FPR(c) = 1 - \text{Specificity}(c) = \frac{FP(c)}{FP(c)+TN(c)} \]

Nilai AUC (area under the curve):

\[ AUC = \int_0^1 TPR(FPR)\,d(FPR). \]

Threshold optimal dipilih dari data training menggunakan indeks Youden:

\[ J = \text{Sensitivity} + \text{Specificity} - 1; \quad c_{\text{optimal}} = \arg\max_c \{J(c)\}. \]

roc_points <- function(actual, prob) {
  thresholds <- c(Inf, sort(unique(prob), decreasing = TRUE), -Inf)
  out <- lapply(thresholds, function(th) {
    pred        <- as.integer(prob >= th)
    tp          <- sum(pred == 1 & actual == 1)
    tn          <- sum(pred == 0 & actual == 0)
    fp          <- sum(pred == 1 & actual == 0)
    fn          <- sum(pred == 0 & actual == 1)
    sensitivity <- safe_div(tp, tp + fn)
    specificity <- safe_div(tn, tn + fp)
    data.frame(
      threshold   = th,
      sensitivity = sensitivity,
      specificity = specificity,
      fpr         = 1 - specificity,
      youden      = sensitivity + specificity - 1
    )
  })
  bind_rows(out)
}

auc_value <- function(roc_df) {
  roc_sorted <- roc_df %>% arrange(fpr, sensitivity)
  sum(
    diff(roc_sorted$fpr) *
      (head(roc_sorted$sensitivity, -1) +
         tail(roc_sorted$sensitivity, -1)) / 2
  )
}

roc_train <- roc_points(train_data$Y_biner, p_train) %>%
  mutate(data = "Training")
roc_test  <- roc_points(test_data$Y_biner,  p_test)  %>%
  mutate(data = "Testing")

auc_train <- auc_value(roc_train)
auc_test  <- auc_value(roc_test)

optimal_train <- roc_train %>%
  filter(is.finite(threshold)) %>%
  arrange(desc(youden), desc(sensitivity)) %>%
  slice(1)

threshold_opt <- optimal_train$threshold[1]

test_at_opt <- roc_points(test_data$Y_biner, p_test) %>%
  filter(is.finite(threshold)) %>%
  slice_min(abs(threshold - threshold_opt), n = 1, with_ties = FALSE) %>%
  mutate(data = "Testing pada threshold optimal")

auc_table <- data.frame(
  Data = c("Training", "Testing"),
  AUC  = round(c(auc_train, auc_test), 3)
)

threshold_table <- optimal_train %>%
  transmute(
    `Threshold optimal` = round(threshold, 3),
    Sensitivity         = round(sensitivity, 3),
    Specificity         = round(specificity, 3),
    `Indeks Youden`     = round(youden, 3)
  )

knitr::kable(auc_table,       caption = "Nilai AUC")

Nilai AUC

Data	AUC
Training	0.712
Testing	0.750

knitr::kable(threshold_table, caption = "Threshold optimal berdasarkan ROC training")

Threshold optimal berdasarkan ROC training

Threshold optimal	Sensitivity	Specificity	Indeks Youden
0.732	0.622	0.777	0.399

Interpretasi output: AUC testing sebesar 0.75 mengukur kemampuan model membedakan responden FAVOR dan OPPOSE pada data baru. Nilai yang lebih besar dari 0,5 berarti model lebih baik daripada tebakan acak. Threshold optimal dari data training adalah 0.732, dipilih karena menghasilkan nilai indeks Youden terbesar sehingga memberi keseimbangan terbaik antara sensitivity dan specificity.

roc_plot <- bind_rows(roc_train, roc_test)

ggplot(roc_plot, aes(x = fpr, y = sensitivity, color = data)) +
  geom_path(linewidth = 1.1) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed",
              color = "#6c757d") +
  geom_point(
    data        = optimal_train,
    aes(x = fpr, y = sensitivity),
    inherit.aes = FALSE,
    color = "#ffb703", fill = "#fb8500",
    shape = 21, size = 4, stroke = 1.2
  ) +
  geom_point(
    data        = test_at_opt,
    aes(x = fpr, y = sensitivity),
    inherit.aes = FALSE,
    color = "#8338ec", fill = "#3a86ff",
    shape = 24, size = 4, stroke = 1.2
  ) +
  coord_equal() +
  scale_color_manual(values = c("Training" = "#0077b6",
                                "Testing"  = "#e76f51")) +
  labs(
    title    = "Kurva ROC Model Regresi Logistik",
    subtitle = paste0(
      "AUC training = ", round(auc_train, 3),
      "; AUC testing = ", round(auc_test, 3),
      "; threshold optimal = ", round(threshold_opt, 3)
    ),
    x     = "False positive rate (1 - specificity)",
    y     = "Sensitivity / true positive rate",
    color = "Data"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Interpretasi output: Kurva ROC yang semakin menjauhi garis diagonal putus-putus menunjukkan performa klasifikasi yang semakin baik. Titik berwarna kuning-oranye menunjukkan posisi threshold optimal pada kurva training, sedangkan titik ungu-biru menunjukkan posisi yang setara pada kurva testing. Jika kedua kurva berdekatan, model cukup stabil pada data baru dan tidak terlalu overfit pada data training.

1.9 Evaluasi dengan Threshold Optimal

Threshold optimal dari ROC training (0.732) kemudian diterapkan pada data testing.

metrics_compare <- bind_rows(
  classification_metrics(test_data$Y_biner, p_test, threshold = 0.5) %>%
    mutate(aturan = "Threshold 0.50"),
  classification_metrics(test_data$Y_biner, p_test,
                         threshold = threshold_opt) %>%
    mutate(aturan = "Threshold optimal ROC")
) %>%
  select(aturan, everything()) %>%
  format_metrics_indonesia() %>%
  bind_cols(
    `Aturan klasifikasi` = c("Threshold 0.50", "Threshold optimal ROC"), .
  ) %>%
  select(`Aturan klasifikasi`, everything())

confusion_opt <- confusion_matrix(test_data$Y_biner, p_test,
                                   threshold = threshold_opt)

knitr::kable(
  metrics_compare %>% mutate(across(where(is.numeric), round, 3)),
  caption = "Perbandingan metrik testing: threshold 0.50 vs threshold optimal"
)

Perbandingan metrik testing: threshold 0.50 vs threshold optimal

Aturan klasifikasi	Threshold	Akurasi	Error rate	Sensitivity	Specificity	Presisi	NPV	F1-score	Balanced accuracy	FPR	FNR
Threshold 0.50	0.500	0.714	0.286	0.681	0.766	0.817	0.610	0.742	0.723	0.234	0.319
Threshold optimal ROC	0.732	0.672	0.328	0.569	0.830	0.837	0.557	0.678	0.700	0.170	0.431

knitr::kable(
  confusion_opt,
  caption = "Confusion matrix testing pada threshold optimal"
)

Confusion matrix testing pada threshold optimal

	Prediksi FAVOR	Prediksi OPPOSE	Sum
Aktual FAVOR	41	31	72
Aktual OPPOSE	8	39	47
Sum	49	70	119

Interpretasi output: Tabel perbandingan menunjukkan dampak perubahan threshold. Jika threshold optimal lebih rendah dari 0.50, model menjadi lebih agresif mendeteksi responden FAVOR; sensitivity meningkat tetapi specificity dapat turun karena sebagian responden OPPOSE ikut terklasifikasi sebagai FAVOR. Pilihan threshold yang tepat bergantung pada konteks: apakah lebih penting mendeteksi semua pendukung hukuman mati (sensitivity tinggi) atau meminimalkan kesalahan klasifikasi responden penolak (specificity tinggi).

test_prob_plot <- test_data %>%
  mutate(peluang_favor = p_test)

ggplot(test_prob_plot, aes(x = peluang_favor, fill = Y_label)) +
  geom_density(alpha = 0.55, color = "white", linewidth = 0.7) +
  geom_vline(
    xintercept = threshold_opt,
    color      = "#fb8500",
    linewidth  = 1.2,
    linetype   = "dashed"
  ) +
  annotate(
    "label",
    x          = threshold_opt,
    y          = Inf,
    label      = paste0("threshold = ", round(threshold_opt, 3)),
    vjust      = 1.4,
    fill       = "#fff3b0",
    color      = "#5f370e",
    label.size = 0
  ) +
  scale_fill_manual(
    values = c("OPPOSE (Menolak)" = "#2a9d8f",
               "FAVOR (Mendukung)" = "#e76f51")
  ) +
  labs(
    title = "Distribusi Peluang Prediksi FAVOR pada Data Testing",
    x     = "Peluang prediksi mendukung hukuman mati (FAVOR)",
    y     = "Kepadatan",
    fill  = "Status aktual"
  ) +
  theme_minimal(base_size = 12)

Interpretasi output: Grafik kepadatan memperlihatkan sebaran peluang prediksi FAVOR untuk responden yang aktualnya FAVOR (oranye) dan OPPOSE (hijau). Semakin terpisah dua kurva, semakin baik model membedakan kedua kelompok. Garis vertikal menunjukkan threshold optimal; observasi di sebelah kanan garis diklasifikasikan sebagai FAVOR, sedangkan di sebelah kiri diklasifikasikan sebagai OPPOSE. Area tumpang tindih antara dua kurva menunjukkan bagian data yang sulit dibedakan oleh model.

1.10 Interpretasi Hasil

1.10.1 Koefisien Model dan Odds Ratio

Berdasarkan hasil estimasi model regresi logistik biner pada data training (n = 472), persamaan logit yang terbentuk adalah:

\[ \text{logit}(\hat{p}_i) = \ln\left(\frac{\hat{p}_i}{1-\hat{p}_i}\right) = 1{,}680 + 0{,}995 \cdot \text{age}_i - 0{,}352 \cdot \mathbf{1}[\text{sex}_i = \text{FEMALE}] - 1{,}600 \cdot \mathbf{1}[\text{polviews}_i = \text{Liberal}]. \]

Interpretasi masing-masing koefisien adalah sebagai berikut.

1. Intercept (\(\hat{\beta}_0\))

Nilai intercept menghasilkan odds ratio sebesar 5,369 (IK 95%: 2,523 – 11,425; p < 0,001). Ini berarti seorang responden laki-laki berusia 0 tahun dengan pandangan Conservative memiliki odds mendukung hukuman mati sebesar 5,369 kali lebih tinggi dibandingkan menolaknya. Nilai ini bersifat ekstrapolasi jauh di luar rentang data (usia minimum = 18 tahun) sehingga tidak memiliki makna substantif langsung, namun diperlukan untuk mengkalibrasi model secara matematis.

2. Usia (age, \(\hat{\beta}_1 = -0{,}005\), OR = 0,995)

Setiap kenaikan usia satu tahun dikaitkan dengan penurunan odds mendukung hukuman mati sebesar faktor 0,995 (IK 95%: 0,984 – 1,007), dengan p-value = 0,426. Karena interval kepercayaan 95% mencakup nilai 1 dan p-value jauh di atas 0,05, usia tidak terbukti sebagai prediktor yang signifikan secara statistik dalam model ini. Meski arah efeknya mengindikasikan sedikit penurunan odds seiring bertambahnya usia, efek tersebut sangat kecil dan tidak dapat dibedakan dari variasi acak pada data.

3. Jenis kelamin (sexFEMALE, \(\hat{\beta}_2\), OR = 0,703)

Responden perempuan memiliki odds mendukung hukuman mati sebesar 0,703 kali odds responden laki-laki (kategori referensi), atau setara dengan penurunan odds sekitar 29,7% (IK 95%: 0,471 – 1,050; p = 0,085). Meskipun arah efeknya konsisten — perempuan cenderung lebih menolak hukuman mati — interval kepercayaan sedikit melewati nilai 1 dan p-value = 0,085 berada di atas ambang 0,05. Dengan demikian, pengaruh jenis kelamin belum signifikan secara statistik pada tingkat signifikansi 5%, meski mendekati batas signifikansi pada tingkat 10%.

4. Pandangan politik (polviewsLiberal, \(\hat{\beta}_3\), OR = 0,202)

Ini adalah prediktor paling dominan dalam model. Responden dengan pandangan Liberal memiliki odds mendukung hukuman mati hanya sebesar 0,202 kali odds responden Conservative (kategori referensi), artinya odds mereka mendukung hukuman mati 79,8% lebih rendah dibandingkan responden Conservative (IK 95%: 0,134 – 0,303; p < 0,001). Interval kepercayaan yang jauh di bawah 1 dan p-value yang sangat kecil memberikan bukti statistik yang sangat kuat bahwa pandangan politik merupakan penentu utama sikap terhadap hukuman mati.

1.10.2 Kebaikan Model (Goodness of Fit)

Deviasi dan penurunan deviance:

Model berhasil menurunkan null deviance dari 631,227 (model tanpa prediktor) menjadi residual deviance 563,852 (model dengan tiga prediktor), menghasilkan penurunan deviance sebesar 67,375 pada 3 derajat bebas. Penurunan ini signifikan secara statistik (\(\chi^2(3) = 67{,}375\), p < 0,001), yang berarti ketiga prediktor secara bersama-sama memberikan kontribusi nyata dalam menjelaskan variasi sikap terhadap hukuman mati.

AUC (Area Under the ROC Curve):

• AUC training = 0,712
• AUC testing = 0,750

AUC testing sebesar 0,750 berada dalam kategori diskriminasi yang cukup baik (acceptable discrimination). Artinya, jika dipilih secara acak satu responden FAVOR dan satu responden OPPOSE, model akan menetapkan peluang prediksi lebih tinggi untuk responden FAVOR sebanyak 75% dari seluruh pasangan yang mungkin. Menariknya, AUC testing (0,750) sedikit lebih tinggi dari AUC training (0,712), yang menunjukkan bahwa model tidak mengalami overfitting dan generalisasinya ke data baru cukup baik.

1.10.3 Evaluasi Klasifikasi

Pada threshold 0,50:

Metrik	Nilai
Akurasi	71,4%
Sensitivity	68,1%
Specificity	76,6%
Presisi	81,7%
F1-score	74,2%
Balanced accuracy	72,3%

Dari 119 observasi testing, model memprediksi dengan benar 49 dari 72 responden FAVOR (sensitivity 68,1%) dan 36 dari 47 responden OPPOSE (specificity 76,6%). Terdapat 23 false negative (responden FAVOR yang diprediksi OPPOSE) dan 11 false positive (responden OPPOSE yang diprediksi FAVOR). Presisi sebesar 81,7% berarti dari semua prediksi FAVOR, 81,7% di antaranya benar-benar FAVOR.

Pada threshold optimal 0,732 (indeks Youden):

Metrik	Nilai
Akurasi	67,2%
Sensitivity	56,9%
Specificity	83,0%
Presisi	83,7%
F1-score	67,8%
Balanced accuracy	70,0%

Dengan menaikkan threshold ke 0,732, model menjadi lebih selektif dalam memprediksi FAVOR. Akibatnya, specificity meningkat dari 76,6% menjadi 83,0% (lebih baik mengenali OPPOSE), namun sensitivity turun dari 68,1% menjadi 56,9% (lebih banyak FAVOR yang terlewat). Untuk konteks penelitian opini publik, threshold 0,50 lebih disarankan karena memberikan keseimbangan yang lebih baik antara sensitivity dan specificity (balanced accuracy 72,3% vs 70,0%) serta F1-score yang lebih tinggi (74,2% vs 67,8%).

1.11 Kesimpulan

Analisis regresi logistik biner terhadap data survei sikap publik mengenai hukuman mati (n = 591) di Amerika Serikat pada tahun 2024 menghasilkan beberapa temuan utama sebagai berikut.

Pertama, mengenai signifikansi prediktor. Dari tiga prediktor yang dimasukkan ke dalam model, hanya pandangan politik (polviews) yang terbukti signifikan secara statistik (p < 0,001). Responden dengan pandangan Liberal memiliki odds mendukung hukuman mati 79,8% lebih rendah dibandingkan responden Conservative (OR = 0,202; IK 95%: 0,134–0,303). Jenis kelamin (sex) menunjukkan arah efek yang konsisten — perempuan cenderung lebih menolak hukuman mati (OR = 0,703) — namun belum mencapai signifikansi statistik pada tingkat 5% (p = 0,085). Usia (age) tidak terbukti berpengaruh signifikan terhadap sikap hukuman mati dalam model ini (OR = 0,995; p = 0,426).

Kedua, mengenai kebaikan model. Model secara keseluruhan signifikan berdasarkan uji penurunan deviance (\(\Delta D = 67{,}375\); df = 3; p < 0,001). Nilai AUC testing sebesar 0,750 menunjukkan kemampuan diskriminasi model termasuk kategori cukup baik, dan tidak terdapat indikasi overfitting karena AUC testing lebih tinggi dari AUC training (0,712).

Ketiga, mengenai performa klasifikasi. Pada threshold 0,50, model mencapai akurasi 71,4% pada data testing, dengan sensitivity 68,1% dan specificity 76,6%. Threshold ini dinilai lebih sesuai dibandingkan threshold optimal Youden (0,732) karena menghasilkan balanced accuracy dan F1-score yang lebih tinggi. Artinya, model lebih mampu mendeteksi responden yang benar-benar mendukung hukuman mati tanpa terlalu banyak salah mengklasifikasikan responden yang menolak.

Secara keseluruhan, pandangan politik merupakan faktor penentu utama sikap terhadap hukuman mati dalam dataset ini. Model regresi logistik biner yang dibangun cukup andal untuk tujuan prediksi awal, meski akurasi yang belum sempurna mencerminkan kompleksitas sikap manusia yang tidak seluruhnya dapat dijelaskan oleh tiga variabel demografis saja. Pengembangan model ke depan dapat mempertimbangkan penambahan variabel seperti tingkat pendidikan, afiliasi agama, atau pengalaman langsung dengan tindak kejahatan untuk meningkatkan daya prediksi model.

2 Regresi Logistik Multinomial

Regresi logistik multinomial adalah model regresi untuk respons kategorik dengan lebih dari dua kategori, ketika kategori tersebut bersifat nominal, yaitu tidak memiliki urutan yang bermakna.

Misalkan respons: \[Y_i \in \{1, 2, \ldots, K\}\] dengan \(K\) kategori. Salah satu kategori dipilih sebagai kategori referensi. Misalkan kategori ke-\(K\) menjadi referensi. Untuk setiap kategori \(k = 1, 2, \ldots, K-1\), model ditulis sebagai:

\[ \log\left(\frac{P(Y_i = k \mid \mathbf{x}_i)}{P(Y_i = K \mid \mathbf{x}_i)}\right) = \beta_{0k} + \beta_{1k} x_{i1} + \cdots + \beta_{pk} x_{ip}. \]

Artinya, model membandingkan peluang masuk kategori \(k\) terhadap peluang masuk kategori referensi. Kasus ini berbeda dari regresi biner: alih-alih satu persamaan logit, model multinomial membangun \((K-1)\) persamaan logit secara simultan, masing-masing membandingkan satu kategori terhadap kategori referensi.

Studi kasus: Analisis faktor-faktor yang memengaruhi status pekerjaan responden berdasarkan data General Social Survey (GSS) 2024.

Tujuan analisis: Membangun model regresi logistik multinomial untuk memprediksi apakah seorang responden bekerja penuh waktu, paruh waktu, mengurus rumah tangga, atau menganggur berdasarkan usia, status pernikahan, tingkat pendidikan, dan jenis kelamin.

Kategori variabel respons (wrkstat):

Kode	Kategori	Keterangan
FT	Working full time	Bekerja penuh waktu
PT	Working part time	Bekerja paruh waktu
HK	Keeping house	Mengurus rumah tangga
UN	Unemployed	Menganggur / sedang mencari kerja

Kategori referensi yang dipilih: Working full time (FT)

2.1 Ringkasan Masalah

Status pekerjaan merupakan variabel multikategori nominal — tidak ada urutan yang inheren antara “bekerja penuh waktu”, “paruh waktu”, “mengurus rumah tangga”, dan “menganggur”. Kasus ini cocok dianalisis dengan regresi logistik multinomial karena variabel respons memiliki lebih dari dua kategori tanpa urutan alami:

• Full Time = Working full time (bekerja penuh waktu)
• Part Time = Working part time (bekerja paruh waktu)
• Keeping House = Keeping house (mengurus rumah tangga)
• Unemployed = Unemployed, laid off, looking for work (menganggur)

Model yang digunakan adalah:

\[ \log\left(\frac{P(Y_i = k \mid X_i)}{P(Y_i = \text{FT} \mid X_i)}\right) = \beta_{0k} + \beta_{1k} X_{1i} + \beta_{2k} X_{2i} + \cdots + \beta_{pk} X_{pi}, \]

dengan \(k \in \{\text{Part Time},\ \text{Keeping House},\ \text{Unemployed}\}\) dan Working full time (FT) sebagai kategori referensi. Model membangun \((K-1) = 3\) persamaan logit secara simultan, masing-masing membandingkan satu kategori terhadap Full Time.

2.2 Sumber Data

Dataset berasal dari General Social Survey (GSS) 2024, survei sosial nasional Amerika Serikat yang diselenggarakan secara berkala oleh NORC di University of Chicago. Dataset terdiri dari 1.575 responden dengan 6 variabel.

2.3 Persiapan Data

2.3.1 Membaca Data dari Excel

raw_data <- readxl::read_excel("Data Multinomial.xlsx", sheet = 1)

ringkasan_data <- data.frame(
  Keterangan = c("Jumlah observasi", "Jumlah variabel"),
  Nilai      = c(nrow(raw_data), ncol(raw_data))
)

knitr::kable(ringkasan_data, caption = "Ukuran dataset GSS 2024")

Ukuran dataset GSS 2024

Keterangan	Nilai
Jumlah observasi	1575
Jumlah variabel	6

Interpretasi output: Dataset berisi 1575 observasi dan 6 variabel. Ukuran sampel ini sangat memadai untuk regresi logistik multinomial dengan empat prediktor.

2.3.2 Penjelasan Variabel

penjelasan <- data.frame(
  `Nama kolom`          = c("id_", "wrkstat", "marital",
                             "age", "degree", "sex"),
  `Keterangan`          = c(
    "Nomor identifikasi responden",
    "Status pekerjaan responden (variabel respons)",
    "Status pernikahan responden",
    "Usia responden (tahun)",
    "Tingkat pendidikan tertinggi",
    "Jenis kelamin responden"
  ),
  `Kategori / Rentang`  = c(
    "—",
    "Working full time, Working part time, Keeping house, Unemployed",
    "Married, Never married",
    "18 – 86 tahun",
    "Less than HS, High school, Associate, Bachelor's, Graduate",
    "MALE, FEMALE"
  ),
  `Tipe dalam analisis` = c(
    "ID (tidak dipakai dalam model)",
    "Respons multinomial (4 kategori)",
    "Kategorik (2 level)",
    "Numerik",
    "Kategorik ordinal (5 level)",
    "Kategorik (2 level)"
  ),
  check.names = FALSE
)

knitr::kable(penjelasan, caption = "Penjelasan variabel dataset GSS 2024")

Penjelasan variabel dataset GSS 2024

Nama kolom	Keterangan	Kategori / Rentang	Tipe dalam analisis
id_	Nomor identifikasi responden	—	ID (tidak dipakai dalam model)
wrkstat	Status pekerjaan responden (variabel respons)	Working full time, Working part time, Keeping house, Unemployed	Respons multinomial (4 kategori)
marital	Status pernikahan responden	Married, Never married	Kategorik (2 level)
age	Usia responden (tahun)	18 – 86 tahun	Numerik
degree	Tingkat pendidikan tertinggi	Less than HS, High school, Associate, Bachelor’s, Graduate	Kategorik ordinal (5 level)
sex	Jenis kelamin responden	MALE, FEMALE	Kategorik (2 level)

Interpretasi output: Terdapat satu variabel respons multinomial (wrkstat) dengan empat kategori, satu prediktor numerik (age), dan tiga prediktor kategorik (marital, degree, sex). Variabel degree memiliki lima level yang diurutkan dari pendidikan terendah ke tertinggi, meski dalam model multinomial tidak diasumsikan ada urutan.

2.3.3 Transformasi dan Pembersihan Data

data_bersih <- raw_data |>
  mutate(
    age = as.numeric(age),

    # Respons: factor dengan referensi "Working full time"
    wrkstat = factor(wrkstat,
                     levels = c("Working full time",
                                "Working part time",
                                "Keeping house",
                                "Unemployed, laid off, looking for work")),

    # Label pendek untuk visualisasi
    wrkstat_label = factor(
      case_when(
        wrkstat == "Working full time"                        ~ "Full Time",
        wrkstat == "Working part time"                        ~ "Part Time",
        wrkstat == "Keeping house"                            ~ "Keeping House",
        wrkstat == "Unemployed, laid off, looking for work"   ~ "Unemployed"
      ),
      levels = c("Full Time", "Part Time", "Keeping House", "Unemployed")
    ),

    # Pendidikan: factor terurut (digunakan sebagai kategorik dalam model)
    degree = factor(degree,
                    levels = c("Less than high school",
                               "High school",
                               "Associate/junior college",
                               "Bachelor's",
                               "Graduate")),

    # Jenis kelamin: MALE sebagai referensi
    sex = factor(sex, levels = c("MALE", "FEMALE")),

    # Status pernikahan: Married sebagai referensi
    marital = factor(marital, levels = c("Married", "Never married"))
  ) |>
  tidyr::drop_na()

knitr::kable(
  head(data_bersih |>
    select(id_, age, sex, marital, degree, wrkstat_label), 8),
  caption = "Contoh data setelah transformasi"
)

Contoh data setelah transformasi

id_	age	sex	marital	degree	wrkstat_label
1	33	MALE	Never married	Bachelor’s	Full Time
4	19	MALE	Never married	High school	Part Time
6	53	MALE	Married	Associate/junior college	Unemployed
7	48	FEMALE	Married	High school	Full Time
9	60	FEMALE	Married	Associate/junior college	Full Time
11	23	FEMALE	Never married	Bachelor’s	Full Time
16	31	FEMALE	Married	Graduate	Full Time
17	64	MALE	Married	Bachelor’s	Full Time

Interpretasi output: Delapan baris pertama data setelah transformasi. Kategori referensi yang dipilih adalah: Working full time untuk wrkstat, MALE untuk sex, Married untuk marital, dan Less than high school untuk degree. Semua perbandingan odds ratio akan mengacu pada kategori-kategori referensi ini.

2.4 Distribusi Kategori Respons

dist_wrkstat <- data_bersih |>
  count(wrkstat_label, name = "Jumlah") |>
  mutate(Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1)) |>
  rename(`Status pekerjaan` = wrkstat_label)

knitr::kable(dist_wrkstat,
             caption = "Distribusi status pekerjaan responden")

Distribusi status pekerjaan responden

Status pekerjaan	Jumlah	Proporsi
Full Time	1045	66.3%
Part Time	229	14.5%
Keeping House	197	12.5%
Unemployed	104	6.6%

Interpretasi output: Mayoritas responden bekerja penuh waktu (Full Time), yang memperkuat alasan pemilihan kategori ini sebagai referensi. Kategori Unemployed merupakan kelompok terkecil, sehingga estimasi parameter untuk kategori ini perlu diinterpretasikan dengan kehati-hatian karena didukung oleh sampel yang lebih kecil.

warna_kategori <- c(
  "Full Time"    = "#2a9d8f",
  "Part Time"    = "#e9c46a",
  "Keeping House"= "#f4a261",
  "Unemployed"   = "#e76f51"
)

ggplot(data_bersih, aes(x = wrkstat_label, fill = wrkstat_label)) +
  geom_bar(width = 0.65, color = "white", linewidth = 0.8) +
  geom_text(stat = "count",
            aes(label = after_stat(count)),
            vjust = -0.4, fontface = "bold") +
  scale_fill_manual(values = warna_kategori) +
  labs(title = "Distribusi Status Pekerjaan Responden GSS 2024",
       x = NULL, y = "Jumlah responden") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

Interpretasi output: Grafik memperlihatkan dominasi responden yang bekerja penuh waktu (Full Time). Ketidakseimbangan ini wajar dalam survei populasi umum dan menjadi alasan kuat menjadikan Full Time sebagai kategori referensi dalam model multinomial.

2.5 Eksplorasi Data

2.5.1 Distribusi Status Pekerjaan Menurut Jenis Kelamin

prop_sex <- data_bersih |>
  count(sex, wrkstat_label) |>
  group_by(sex) |>
  mutate(prop = n / sum(n),
         label = paste0(round(prop * 100, 1), "%")) |>
  ungroup()

ggplot(prop_sex, aes(x = sex, y = prop, fill = wrkstat_label)) +
  geom_col(width = 0.65, color = "white", linewidth = 0.6) +
  geom_text(aes(label = label),
            position = position_stack(vjust = 0.5),
            color = "white", fontface = "bold", size = 3.2) +
  scale_fill_manual(values = warna_kategori) +
  scale_y_continuous(labels = scales::percent) +
  labs(title    = "Proporsi Status Pekerjaan Menurut Jenis Kelamin",
       x = "Jenis kelamin", y = "Proporsi", fill = "Status pekerjaan") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Interpretasi output: Grafik stacked bar menampilkan perbedaan komposisi status pekerjaan antara laki-laki dan perempuan. Jika proporsi Keeping House jauh lebih tinggi pada perempuan, ini mengindikasikan bahwa jenis kelamin merupakan prediktor penting khususnya untuk kategori tersebut.

2.5.2 Distribusi Status Pekerjaan Menurut Tingkat Pendidikan

prop_degree <- data_bersih |>
  count(degree, wrkstat_label) |>
  group_by(degree) |>
  mutate(prop = n / sum(n),
         label = ifelse(prop > 0.05,
                        paste0(round(prop * 100, 0), "%"), "")) |>
  ungroup()

ggplot(prop_degree, aes(x = degree, y = prop, fill = wrkstat_label)) +
  geom_col(width = 0.72, color = "white", linewidth = 0.6) +
  geom_text(aes(label = label),
            position = position_stack(vjust = 0.5),
            color = "white", fontface = "bold", size = 3.0) +
  scale_fill_manual(values = warna_kategori) +
  scale_y_continuous(labels = scales::percent) +
  scale_x_discrete(labels = c(
    "Less than high school"    = "< SMA",
    "High school"              = "SMA",
    "Associate/junior college" = "D3",
    "Bachelor's"               = "S1",
    "Graduate"                 = "S2/S3"
  )) +
  labs(title    = "Proporsi Status Pekerjaan Menurut Tingkat Pendidikan",
       x = "Tingkat pendidikan", y = "Proporsi",
       fill = "Status pekerjaan") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Interpretasi output: Grafik ini menampilkan bagaimana komposisi status pekerjaan berubah seiring meningkatnya tingkat pendidikan. Jika Full Time meningkat dan Unemployed menurun seiring kenaikan pendidikan, ini menjadi indikasi kuat bahwa degree merupakan prediktor relevan dalam model.

2.5.3 Distribusi Usia Menurut Status Pekerjaan

ggplot(data_bersih, aes(x = wrkstat_label, y = age, fill = wrkstat_label)) +
  geom_boxplot(alpha = 0.75, color = "#243142",
               outlier.alpha = 0.4, outlier.size = 1.2) +
  scale_fill_manual(values = warna_kategori) +
  labs(title = "Distribusi Usia Menurut Status Pekerjaan",
       x = NULL, y = "Usia (tahun)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

Interpretasi output: Boxplot memperlihatkan sebaran usia untuk setiap kategori status pekerjaan. Perbedaan median usia yang mencolok antar kategori mengindikasikan bahwa age berperan dalam menentukan status pekerjaan seseorang.

2.5.4 Ringkasan Numerik Usia

usia_summary <- data_bersih |>
  group_by(wrkstat_label) |>
  summarise(
    Jumlah   = n(),
    Rata2    = mean(age, na.rm = TRUE),
    Median   = median(age, na.rm = TRUE),
    Std      = sd(age, na.rm = TRUE),
    Min      = min(age),
    Maks     = max(age),
    .groups  = "drop"
  ) |>
  rename(`Status pekerjaan` = wrkstat_label) |>
  mutate(across(where(is.numeric), round, 2))

knitr::kable(usia_summary,
             caption = "Ringkasan usia menurut status pekerjaan")

Ringkasan usia menurut status pekerjaan

Status pekerjaan	Jumlah	Rata2	Median	Std	Min	Maks
Full Time	1045	41.39	40.0	12.56	18	85
Part Time	229	45.77	44.0	17.38	18	85
Keeping House	197	42.45	39.0	14.73	18	86
Unemployed	104	38.27	36.5	14.00	18	73

Interpretasi output: Tabel merangkum statistik deskriptif usia untuk setiap kategori status pekerjaan. Perbedaan rata-rata dan median usia antar kategori memberikan gambaran awal mengenai peran usia sebagai prediktor.

2.6 Estimasi Model Multinomial dengan R

Model diestimasi menggunakan fungsi multinom() dari paket nnet. Kategori referensi untuk wrkstat adalah Working full time.

model_multi <- nnet::multinom(
  wrkstat ~ age + sex + marital + degree,
  data  = data_bersih,
  trace = FALSE      # sembunyikan iterasi konvergensi
)

summary(model_multi)

## Call:
## nnet::multinom(formula = wrkstat ~ age + sex + marital + degree, 
##     data = data_bersih, trace = FALSE)
## 
## Coefficients:
##                                        (Intercept)         age sexFEMALE
## Working part time                        -2.901070 0.030762276 0.8249242
## Keeping house                            -1.792936 0.005173082 2.2162268
## Unemployed, laid off, looking for work   -2.267533 0.005188112 0.4806244
##                                        maritalNever married degreeHigh school
## Working part time                                 0.2898141        -0.1862848
## Keeping house                                    -0.6633728        -0.8380199
## Unemployed, laid off, looking for work            1.4435573        -1.2044847
##                                        degreeAssociate/junior college
## Working part time                                          -0.5425302
## Keeping house                                              -2.2368218
## Unemployed, laid off, looking for work                     -1.6517380
##                                        degreeBachelor's degreeGraduate
## Working part time                            -0.8297569      -1.312131
## Keeping house                                -1.8456750      -2.499093
## Unemployed, laid off, looking for work       -2.0342674      -2.563259
## 
## Std. Errors:
##                                        (Intercept)         age sexFEMALE
## Working part time                        0.4195370 0.005717080 0.1539669
## Keeping house                            0.4368631 0.006603117 0.2174324
## Unemployed, laid off, looking for work   0.5388935 0.008461851 0.2159800
##                                        maritalNever married degreeHigh school
## Working part time                                 0.1653295         0.2855703
## Keeping house                                     0.1888033         0.2671407
## Unemployed, laid off, looking for work            0.2766001         0.2927154
##                                        degreeAssociate/junior college
## Working part time                                           0.3544072
## Keeping house                                               0.4321605
## Unemployed, laid off, looking for work                      0.4484440
##                                        degreeBachelor's degreeGraduate
## Working part time                             0.3096431      0.3525457
## Keeping house                                 0.3067165      0.3712005
## Unemployed, laid off, looking for work        0.3732728      0.5221890
## 
## Residual Deviance: 2775.321 
## AIC: 2823.321

Interpretasi output: multinom() tidak langsung memberikan p-value. Koefisien pada baris Part Time, Keeping House, dan Unemployed masing-masing menjelaskan perubahan log-odds masuk kategori tersebut dibandingkan Full Time.

2.7 Ringkasan Koefisien, Standard Error, z-value, dan p-value

Fungsi multinom() tidak langsung memberikan p-value. Kita dapat menghitungnya dari:

\[z = \frac{\hat{\beta}}{SE(\hat{\beta})}, \qquad p\text{-value} = 2\{1 - \Phi(|z|)\}.\]

# Ekstrak koefisien, SE, z, p-value
coef_mat  <- summary(model_multi)$coefficients
se_mat    <- summary(model_multi)$standard.errors
z_mat     <- coef_mat / se_mat
p_mat     <- 2 * (1 - pnorm(abs(z_mat)))

# Susun ke data frame panjang
tidy_multi <- function(coef, se, z, p) {
  kategori  <- rownames(coef)
  prediktor <- colnames(coef)
  out <- expand.grid(Kategori = kategori, Prediktor = prediktor,
                     stringsAsFactors = FALSE)
  out$Koefisien <- as.vector(t(coef))
  out$SE        <- as.vector(t(se))
  out$z         <- as.vector(t(z))
  out$p_value   <- as.vector(t(p))
  out$RRR       <- exp(out$Koefisien)
  out$CI_low    <- exp(out$Koefisien - 1.96 * out$SE)
  out$CI_high   <- exp(out$Koefisien + 1.96 * out$SE)
  out
}

rrr_df <- tidy_multi(coef_mat, se_mat, z_mat, p_mat) |>
  filter(Prediktor != "(Intercept)") |>
  mutate(
    Signifikan = case_when(
      p_value < 0.001 ~ "***",
      p_value < 0.01  ~ "**",
      p_value < 0.05  ~ "*",
      p_value < 0.1   ~ ".",
      TRUE            ~ ""
    )
  )

# Label pendek untuk kategori
label_kat <- c(
  "Working part time"                      = "Part Time vs Full Time",
  "Keeping house"                          = "Keeping House vs Full Time",
  "Unemployed, laid off, looking for work" = "Unemployed vs Full Time"
)

rrr_display <- rrr_df |>
  mutate(
    Perbandingan = label_kat[Kategori],
    `RRR`        = round(RRR, 4),
    `SE`         = round(SE, 4),
    `z-value`    = round(z, 4),
    `p-value`    = signif(p_value, 3),
    `IK 95%`     = paste0(round(CI_low, 3), " – ", round(CI_high, 3))
  ) |>
  select(Perbandingan, Prediktor, `RRR`, `SE`, `z-value`, `p-value`,
         `IK 95%`, Signifikan)

knitr::kable(rrr_display,
             caption = "Ringkasan hasil regresi logistik multinomial (referensi: Full Time)")

Ringkasan hasil regresi logistik multinomial (referensi: Full Time)

Perbandingan	Prediktor	RRR	SE	z-value	p-value	IK 95%	Signifikan
Part Time vs Full Time	age	1.3362	0.1653	1.7529	7.96e-02	0.966 – 1.848	.
Keeping House vs Full Time	age	0.8300	0.2856	-0.6523	5.14e-01	0.474 – 1.453
Unemployed vs Full Time	age	0.5813	0.3544	-1.5308	1.26e-01	0.29 – 1.164
Part Time vs Full Time	sexFEMALE	0.4362	0.3096	-2.6797	7.37e-03	0.238 – 0.8	**
Keeping House vs Full Time	sexFEMALE	0.2692	0.3525	-3.7219	1.98e-04	0.135 – 0.537	***
Unemployed vs Full Time	sexFEMALE	0.1665	0.4369	-4.1041	4.06e-05	0.071 – 0.392	***
Part Time vs Full Time	maritalNever married	1.0052	0.0066	0.7834	4.33e-01	0.992 – 1.018
Keeping House vs Full Time	maritalNever married	9.1727	0.2174	10.1927	0.00e+00	5.99 – 14.047	***
Unemployed vs Full Time	maritalNever married	0.5151	0.1888	-3.5136	4.42e-04	0.356 – 0.746	***
Part Time vs Full Time	degreeHigh school	0.4326	0.2671	-3.1370	1.71e-03	0.256 – 0.73	**
Keeping House vs Full Time	degreeHigh school	0.1068	0.4322	-5.1759	2.00e-07	0.046 – 0.249	***
Unemployed vs Full Time	degreeHigh school	0.1579	0.3067	-6.0175	0.00e+00	0.087 – 0.288	***
Part Time vs Full Time	degreeAssociate/junior college	0.0822	0.3712	-6.7325	0.00e+00	0.04 – 0.17	***
Keeping House vs Full Time	degreeAssociate/junior college	0.1036	0.5389	-4.2078	2.58e-05	0.036 – 0.298	***
Unemployed vs Full Time	degreeAssociate/junior college	1.0052	0.0085	0.6131	5.40e-01	0.989 – 1.022
Part Time vs Full Time	degreeBachelor’s	1.6171	0.2160	2.2253	2.61e-02	1.059 – 2.469	*
Keeping House vs Full Time	degreeBachelor’s	4.2357	0.2766	5.2189	2.00e-07	2.463 – 7.284	***
Unemployed vs Full Time	degreeBachelor’s	0.2998	0.2927	-4.1149	3.87e-05	0.169 – 0.532	***
Part Time vs Full Time	degreeGraduate	0.1917	0.4484	-3.6833	2.30e-04	0.08 – 0.462	***
Keeping House vs Full Time	degreeGraduate	0.1308	0.3733	-5.4498	1.00e-07	0.063 – 0.272	***
Unemployed vs Full Time	degreeGraduate	0.0771	0.5222	-4.9087	9.00e-07	0.028 – 0.214	***

Catatan: *** p < 0,001 | ** p < 0,01 | * p < 0,05 | . p < 0,1

Interpretasi output: Tabel menampilkan tiga blok perbandingan. Setiap baris adalah satu koefisien prediktor untuk satu perbandingan kategori.

2.8 Interpretasi Koefisien Multinomial

Misalkan kategori referensi adalah Working full time. Maka:

• Koefisien pada baris Part Time menjelaskan perubahan log-odds bekerja paruh waktu dibanding bekerja penuh waktu.
• Koefisien pada baris Keeping House menjelaskan perubahan log-odds mengurus rumah tangga dibanding bekerja penuh waktu.
• Koefisien pada baris Unemployed menjelaskan perubahan log-odds menganggur dibanding bekerja penuh waktu.

Jika \(\exp(\beta) > 1\), maka kenaikan prediktor meningkatkan odds relatif masuk kategori tersebut dibanding kategori referensi. Jika \(\exp(\beta) < 1\), maka sebaliknya.

Contoh interpretasi: Jika koefisien sexFEMALE untuk kategori Keeping House bernilai positif dan besar, maka responden perempuan memiliki kecenderungan yang secara substansial lebih besar untuk berada dalam status mengurus rumah tangga dibandingkan laki-laki, dengan asumsi variabel lain konstan.

2.8.1 Visualisasi RRR

rrr_plot <- rrr_df |>
  filter(!grepl("Intercept", Prediktor)) |>
  mutate(
    Perbandingan = label_kat[Kategori],
    Prediktor    = gsub("degree|sex|marital", "", Prediktor)
  )

ggplot(rrr_plot,
       aes(x = RRR, y = reorder(Prediktor, RRR),
           color = Perbandingan)) +
  geom_vline(xintercept = 1, linetype = "dashed",
             color = "#6c757d", linewidth = 0.8) +
  geom_errorbarh(aes(xmin = CI_low, xmax = CI_high),
                 height = 0.25, linewidth = 0.7,
                 position = position_dodge(width = 0.6)) +
  geom_point(size = 2.8,
             position = position_dodge(width = 0.6)) +
  scale_color_manual(values = c(
    "Part Time vs Full Time"    = "#e9c46a",
    "Keeping House vs Full Time"= "#f4a261",
    "Unemployed vs Full Time"   = "#e76f51"
  )) +
  scale_x_log10() +
  labs(title    = "Relative Risk Ratio (RRR) Model Multinomial",
       subtitle = "Skala logaritmik; garis putus-putus = RRR 1 (tidak ada efek)",
       x        = "RRR (skala log)",
       y        = NULL,
       color    = "Perbandingan") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Interpretasi output: Forest plot RRR menampilkan estimasi titik dan interval kepercayaan 95% untuk setiap koefisien. Titik di sebelah kanan garis putus-putus (RRR > 1) berarti prediktor meningkatkan kecenderungan masuk kategori tersebut dibanding Full Time. Jika interval kepercayaan tidak melewati garis 1, koefisien signifikan secara statistik.

2.9 Prediksi Probabilitas Multinomial

prob_all <- predict(model_multi, newdata = data_bersih, type = "probs")
pred_kelas <- predict(model_multi, newdata = data_bersih, type = "class")

head(prob_all)

##   Working full time Working part time Keeping house
## 1         0.8528441        0.07539336    0.01369875
## 2         0.7719299        0.08442271    0.03159040
## 3         0.8246349        0.13452384    0.01928554
## 4         0.4225037        0.19255022    0.35772929
## 5         0.5787316        0.26717386    0.12872678
## 6         0.7180141        0.10647744    0.10045527
##   Unemployed, laid off, looking for work
## 1                             0.05806376
## 2                             0.11205702
## 3                             0.02155570
## 4                             0.02721680
## 5                             0.02536777
## 6                             0.07505322

Tambahkan prediksi ke data:

data_bersih_pred <- data_bersih |>
  dplyr::bind_cols(as.data.frame(prob_all)) |>
  mutate(
    prediksi = factor(
      case_when(
        pred_kelas == "Working full time"                      ~ "Full Time",
        pred_kelas == "Working part time"                      ~ "Part Time",
        pred_kelas == "Keeping house"                          ~ "Keeping House",
        pred_kelas == "Unemployed, laid off, looking for work" ~ "Unemployed"
      ),
      levels = c("Full Time", "Part Time", "Keeping House", "Unemployed")
    )
  )

# Tampilkan contoh prediksi
pred_preview <- head(
  data.frame(
    `Status aktual`  = data_bersih_pred$wrkstat_label,
    `Pred. kelas`    = data_bersih_pred$prediksi,
    `P(Full Time)`   = round(prob_all[, "Working full time"], 3),
    `P(Part Time)`   = round(prob_all[, "Working part time"], 3),
    `P(Keep. House)` = round(prob_all[, "Keeping house"], 3),
    `P(Unemployed)`  = round(prob_all[,
      "Unemployed, laid off, looking for work"], 3),
    check.names = FALSE
  ),
  8
)

knitr::kable(pred_preview,
             caption = "Contoh prediksi probabilitas (6 observasi pertama)")

Contoh prediksi probabilitas (6 observasi pertama)

Status aktual	Pred. kelas	P(Full Time)	P(Part Time)	P(Keep. House)	P(Unemployed)
Full Time	Full Time	0.853	0.075	0.014	0.058
Part Time	Full Time	0.772	0.084	0.032	0.112
Unemployed	Full Time	0.825	0.135	0.019	0.022
Full Time	Full Time	0.423	0.193	0.358	0.027
Full Time	Full Time	0.579	0.267	0.129	0.025
Full Time	Full Time	0.718	0.106	0.100	0.075
Full Time	Full Time	0.800	0.070	0.118	0.012
Full Time	Full Time	0.815	0.140	0.030	0.015

Interpretasi output: Setiap baris menampilkan empat peluang prediksi yang jumlahnya selalu = 1. Model mengklasifikasikan observasi ke kategori dengan peluang tertinggi.

2.10 Confusion Matrix Sederhana

conf_multi <- table(
  Aktual   = data_bersih_pred$wrkstat_label,
  Prediksi = data_bersih_pred$prediksi
)

conf_multi

##                Prediksi
## Aktual          Full Time Part Time Keeping House Unemployed
##   Full Time          1035         2             8          0
##   Part Time           219         4             6          0
##   Keeping House       178         6            13          0
##   Unemployed           96         2             6          0

knitr::kable(addmargins(conf_multi),
             caption = "Confusion matrix (baris = aktual, kolom = prediksi)")

Confusion matrix (baris = aktual, kolom = prediksi)

	Full Time	Part Time	Keeping House	Unemployed	Sum
Full Time	1035	2	8	0	1045
Part Time	219	4	6	0	229
Keeping House	178	6	13	0	197
Unemployed	96	2	6	0	104
Sum	1528	14	33	0	1575

accuracy_multi <- sum(diag(conf_multi)) / sum(conf_multi)
accuracy_multi

## [1] 0.6679365

Interpretasi output: Diagonal confusion matrix menunjukkan prediksi yang benar untuk setiap kategori. Akurasi tidak selalu cukup untuk menilai model, khususnya jika kategori tidak seimbang. Untuk analisis yang lebih lengkap dapat digunakan ukuran lain seperti macro-F1, balanced accuracy, atau log-loss.

2.10.1 Metrik Evaluasi Per Kategori

actual_vec <- data_bersih_pred$wrkstat_label
pred_vec   <- data_bersih_pred$prediksi

kelas_levels <- levels(actual_vec)

per_kelas <- lapply(kelas_levels, function(k) {
  tp <- sum(actual_vec == k & pred_vec == k)
  fp <- sum(actual_vec != k & pred_vec == k)
  fn <- sum(actual_vec == k & pred_vec != k)
  tn <- sum(actual_vec != k & pred_vec != k)
  sensitivity <- ifelse((tp + fn) == 0, NA, tp / (tp + fn))
  precision   <- ifelse((tp + fp) == 0, NA, tp / (tp + fp))
  specificity <- ifelse((tn + fp) == 0, NA, tn / (tn + fp))
  f1          <- ifelse(is.na(precision) | is.na(sensitivity) |
                        (precision + sensitivity) == 0, NA,
                        2 * precision * sensitivity /
                          (precision + sensitivity))
  data.frame(Kategori = k,
             TP = tp, FP = fp, FN = fn, TN = tn,
             Sensitivity = round(sensitivity, 3),
             Specificity = round(specificity, 3),
             Precision   = round(precision, 3),
             F1          = round(f1, 3))
}) |> dplyr::bind_rows()

knitr::kable(per_kelas,
             caption = "Metrik evaluasi per kategori")

Metrik evaluasi per kategori

Kategori	TP	FP	FN	TN	Sensitivity	Specificity	Precision	F1
Full Time	1035	493	10	37	0.990	0.070	0.677	0.805
Part Time	4	10	225	1336	0.017	0.993	0.286	0.033
Keeping House	13	20	184	1358	0.066	0.985	0.394	0.113
Unemployed	0	0	104	1471	0.000	1.000	NA	NA

overall_table <- data.frame(
  Metrik = c("Overall Accuracy",
             "Macro-avg Sensitivity",
             "Macro-avg Precision",
             "Macro-avg F1-score"),
  Nilai  = round(c(
    accuracy_multi,
    mean(per_kelas$Sensitivity, na.rm = TRUE),
    mean(per_kelas$Precision,   na.rm = TRUE),
    mean(per_kelas$F1,          na.rm = TRUE)
  ), 3)
)

knitr::kable(overall_table,
             caption = "Metrik evaluasi keseluruhan (macro-average)")

Metrik evaluasi keseluruhan (macro-average)

Metrik	Nilai
Overall Accuracy	0.668
Macro-avg Sensitivity	0.268
Macro-avg Precision	0.452
Macro-avg F1-score	0.317

Interpretasi output — metrik per kategori: Sensitivity mengukur seberapa baik model mengenali observasi yang benar-benar masuk kategori tersebut. Precision mengukur seberapa andal prediksi untuk kategori itu. F1-score menyeimbangkan keduanya.

Interpretasi output — metrik keseluruhan: Macro-average menghitung rata-rata sederhana metrik per kategori tanpa mempertimbangkan ukuran kelas, sehingga memberikan gambaran yang lebih adil terhadap kategori minoritas seperti Unemployed.

2.11 Visualisasi Probabilitas Prediksi

Kita dapat melihat bagaimana probabilitas pilihan status pekerjaan berubah terhadap usia, dengan variabel lain ditahan pada nilai/kategori referensinya.

grid_multi <- expand.grid(
  age     = seq(min(data_bersih$age), max(data_bersih$age), length.out = 100),
  sex     = factor("MALE",        levels = levels(data_bersih$sex)),
  marital = factor("Married",     levels = levels(data_bersih$marital)),
  degree  = factor("High school", levels = levels(data_bersih$degree))
)

grid_prob <- predict(model_multi, newdata = grid_multi, type = "probs")

grid_plot <- grid_multi |>
  dplyr::bind_cols(as.data.frame(grid_prob)) |>
  tidyr::pivot_longer(
    cols      = c("Working full time", "Working part time",
                  "Keeping house",
                  "Unemployed, laid off, looking for work"),
    names_to  = "status",
    values_to = "probabilitas"
  ) |>
  mutate(status = case_when(
    status == "Working full time"                      ~ "Full Time",
    status == "Working part time"                      ~ "Part Time",
    status == "Keeping house"                          ~ "Keeping House",
    status == "Unemployed, laid off, looking for work" ~ "Unemployed"
  ),
  status = factor(status,
                  levels = c("Full Time", "Part Time",
                             "Keeping House", "Unemployed")))

ggplot(grid_plot, aes(x = age, y = probabilitas, color = status)) +
  geom_line(linewidth = 1.35) +
  scale_y_continuous(labels = scales::percent) +
  scale_color_manual(values = warna_kategori) +
  labs(title    = "Prediksi Probabilitas Status Pekerjaan",
       subtitle = "Variabel lain: MALE, Married, High school.",
       x        = "Usia (tahun)",
       y        = "Probabilitas prediksi",
       color    = "Status pekerjaan") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Interpretasi output: Grafik ini membantu membaca hasil model secara lebih substantif. Jika garis Keeping House meningkat seiring bertambahnya usia, model menunjukkan bahwa responden yang lebih tua lebih cenderung berada dalam status mengurus rumah tangga. Interpretasi seperti ini sering lebih mudah dipahami dibanding hanya melihat tabel koefisien.

2.11.1 Visualisasi Confusion Matrix

cm_df <- as.data.frame(table(
  Aktual   = data_bersih_pred$wrkstat_label,
  Prediksi = data_bersih_pred$prediksi
)) |>
  group_by(Aktual) |>
  mutate(prop = Freq / sum(Freq)) |>
  ungroup()

ggplot(cm_df, aes(x = Prediksi, y = Aktual, fill = prop)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(aes(label = paste0(Freq, "\n(", scales::percent(prop, accuracy = 1), ")")),
            fontface = "bold", size = 3.4) +
  scale_fill_gradient(low = "#f8f9fa", high = "#2a9d8f",
                      labels = scales::percent) +
  labs(title    = "Confusion Matrix (Proporsi per Baris)",
       subtitle = "Angka = jumlah observasi; persentase = dari total aktual per kategori",
       x = "Prediksi", y = "Aktual",
       fill = "Proporsi") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

Interpretasi output: Heatmap confusion matrix memudahkan identifikasi pola kesalahan. Sel diagonal yang berwarna lebih gelap menunjukkan performa yang lebih baik untuk kategori tersebut.

2.12 Catatan Metodologis

• Model ini menggunakan baseline-category logit dengan Full Time sebagai referensi. Pilihan referensi tidak mengubah prediksi akhir, tetapi mengubah arah interpretasi RRR.
• Variabel degree diperlakukan sebagai kategorik nominal meskipun memiliki urutan alami. Jika urutan ingin dieksploitasi, dapat dipertimbangkan menggunakan model ordinal atau mengkode degree sebagai skor numerik.
• Akurasi tidak selalu cukup untuk menilai model pada kelas yang tidak seimbang. Ukuran seperti macro-F1, balanced accuracy, atau log-loss memberikan evaluasi yang lebih komprehensif.
• Untuk pengembangan lanjutan, dapat dipertimbangkan penambahan variabel seperti wilayah tempat tinggal, jumlah tanggungan, atau sektor industri untuk meningkatkan akurasi model.

2.13 Kesimpulan

Analisis regresi logistik multinomial terhadap data GSS 2024 menghasilkan beberapa temuan utama:

Pertama, formulasi model. Model membentuk tiga persamaan logit secara simultan — masing-masing membandingkan Part Time, Keeping House, dan Unemployed terhadap Full Time sebagai kategori referensi. Koefisien diinterpretasikan sebagai relative risk ratio (RRR) setelah dieksponensialkan.

Kedua, prediktor dominan. Jenis kelamin (sex) dan tingkat pendidikan (degree) kemungkinan besar merupakan prediktor terkuat, khususnya untuk memisahkan kategori Keeping House dan Full Time, serta Unemployed dan Full Time. Perempuan memiliki kecenderungan lebih tinggi untuk berada dalam status Keeping House dibandingkan laki-laki. Pendidikan yang lebih tinggi umumnya dikaitkan dengan kecenderungan lebih rendah untuk pengangguran.

Ketiga, performa klasifikasi. Overall accuracy pada data sebesar 66.8%. Ketimpangan performa antar kategori mencerminkan ketidakseimbangan jumlah sampel, di mana Full Time mendominasi dataset.

Keempat, keterbatasan. Model ini memberikan gambaran awal yang informatif mengenai faktor demografis yang berkaitan dengan status pekerjaan, namun tidak boleh diinterpretasikan sebagai hubungan kausal. Asumsi IIA yang melekat pada model multinomial perlu diverifikasi secara teoritis sesuai konteks.

3 Regresi Logistik Ordinal

Regresi logistik ordinal digunakan ketika variabel respons bersifat kategorik dengan urutan alami. Berbeda dari regresi multinomial yang memperlakukan semua kategori setara, model ordinal mengeksploitasi informasi urutan tersebut sehingga lebih efisien secara statistik dan lebih kaya secara interpretatif.

Kasus ini berbeda dari regresi biner maupun multinomial:

Model	Tipe respons	Contoh
Logistik biner	2 kategori tanpa urutan	Setuju / Tidak setuju
Logistik multinomial	≥ 3 kategori tanpa urutan	Jenis pekerjaan
Logistik ordinal	≥ 3 kategori dengan urutan	Tingkat kebahagiaan

Studi kasus: Analisis determinan tingkat kebahagiaan responden berdasarkan data General Social Survey (GSS) 2024.

Tujuan analisis: Membangun model regresi logistik ordinal untuk memprediksi tingkat kebahagiaan (Not too happy < Pretty happy < Very happy) berdasarkan status pernikahan, kondisi kesehatan, dan pendapatan riil responden.

3.1 Ringkasan Masalah

Tingkat kebahagiaan subjektif merupakan variabel ordinal — kategorinya memiliki urutan yang jelas namun jarak antar-kategori tidak diasumsikan sama. Model yang tepat adalah regresi logistik ordinal, bukan multinomial (yang mengabaikan urutan) maupun regresi linear biasa (yang mengasumsikan jarak sama antar-kategori).

Model yang digunakan adalah:

\[ \log\left( \frac{P(Y_i \leq j \mid \mathbf{x}_i)}{P(Y_i > j \mid \mathbf{x}_i)} \right) = \alpha_j + \beta_1 \cdot \text{log\_income}_i + \beta_2 \cdot \mathbf{1}[\text{marital}_i = \text{Married}] + \beta_3 \cdot \mathbf{1}[\text{health}_i = \text{Fair}] + \cdots + \beta_5 \cdot \mathbf{1}[\text{health}_i = \text{Excellent}], \]

dengan \(j \in \{1, 2\}\) (dua batas kumulatif untuk tiga kategori) dan vektor koefisien \(\boldsymbol{\beta}\) yang sama untuk semua \(j\) (proportional odds assumption). Model membangun \(J - 1 = 2\) persamaan cumulative logit secara simultan, hanya dibedakan oleh cutpoint \(\alpha_1 < \alpha_2\).

3.2 Sumber Data

Dataset berasal dari General Social Survey (GSS) 2024, survei sosial nasional Amerika Serikat yang diselenggarakan oleh NORC di University of Chicago. Dataset terdiri dari 2.124 responden dengan 5 variabel.

3.3 Persiapan Data

3.3.1 Membaca Data dari Excel

raw_data <- readxl::read_excel("Data Ordinal.xlsx", sheet = 1)

data.frame(
  Keterangan = c("Jumlah observasi", "Jumlah variabel"),
  Nilai      = c(nrow(raw_data), ncol(raw_data))
) |>
  kable(caption = "Ukuran dataset GSS 2024") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE)

Ukuran dataset GSS 2024

Keterangan	Nilai
Jumlah observasi	2124
Jumlah variabel	5

Interpretasi output: Dataset berisi 2124 observasi dan 5 variabel. Ukuran sampel ini memadai untuk regresi logistik ordinal dengan tiga prediktor.

3.3.2 Penjelasan Variabel

kamus <- data.frame(
  `Nama kolom`          = c("id_","marital","happy","health","realinc"),
  `Keterangan`          = c(
    "Nomor identifikasi responden",
    "Status pernikahan responden",
    "Tingkat kebahagiaan responden (variabel respons)",
    "Kondisi kesehatan responden",
    "Pendapatan riil keluarga (USD per tahun)"
  ),
  `Kategori / Rentang`  = c(
    "—",
    "Married, Never married",
    "Not too happy < Pretty happy < Very happy",
    "Poor < Fair < Good < Excellent",
    "Kontinu; min ~181, maks ~117.565"
  ),
  `Tipe dalam analisis` = c(
    "ID (tidak dipakai dalam model)",
    "Kategorik (2 level)",
    "Respons ordinal (3 kategori)",
    "Kategorik ordinal (4 level, dipakai sebagai kategorik)",
    "Numerik (log-transformasi disarankan)"
  ),
  check.names = FALSE
)

kamus |>
  kable(caption = "Penjelasan variabel dataset GSS 2024") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = TRUE)

Penjelasan variabel dataset GSS 2024

Nama kolom	Keterangan	Kategori / Rentang	Tipe dalam analisis
id_	Nomor identifikasi responden	—	ID (tidak dipakai dalam model)
marital	Status pernikahan responden	Married, Never married	Kategorik (2 level)
happy	Tingkat kebahagiaan responden (variabel respons)	Not too happy < Pretty happy < Very happy	Respons ordinal (3 kategori)
health	Kondisi kesehatan responden	Poor < Fair < Good < Excellent	Kategorik ordinal (4 level, dipakai sebagai kategorik)
realinc	Pendapatan riil keluarga (USD per tahun)	Kontinu; min ~181, maks ~117.565	Numerik (log-transformasi disarankan)

Interpretasi output: Variabel respons happy memiliki tiga kategori berurutan sehingga tepat dimodelkan dengan regresi logistik ordinal. Variabel health juga bersifat ordinal, tetapi diperlakukan sebagai prediktor kategorik dalam model ini. Pendapatan (realinc) memiliki distribusi yang sangat miring ke kanan sehingga perlu ditransformasi logaritmik sebelum dimasukkan ke model.

3.3.3 Transformasi dan Pembersihan Data

data_bersih <- raw_data |>
  mutate(
    realinc = as.numeric(realinc),

    # Log-transformasi pendapatan untuk menstabilkan distribusi
    log_income = log(realinc + 1),

    # Respons ordinal: urutan dari rendah ke tinggi
    happy = ordered(happy,
                    levels = c("Not too happy", "Pretty happy", "Very happy")),

    # Kondisi kesehatan sebagai factor urutan; Poor sebagai referensi
    health = factor(health,
                    levels = c("Poor","Fair","Good","Excellent")),

    # Status pernikahan; Never married sebagai referensi
    marital = factor(marital,
                     levels = c("Never married", "Married"))
  ) |>
  tidyr::drop_na()

head(data_bersih |>
  select(id_, marital, health, realinc, log_income, happy), 8) |>
  kable(caption = "Contoh data setelah transformasi",
        digits = 2) |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE)

Contoh data setelah transformasi

id_	marital	health	realinc	log_income	happy
1	Never married	Good	43560.0	10.68	Not too happy
2	Never married	Good	24502.5	10.11	Pretty happy
3	Married	Good	43560.0	10.68	Pretty happy
4	Never married	Good	58080.0	10.97	Pretty happy
6	Married	Good	2722.5	7.91	Pretty happy
7	Married	Good	117564.8	11.67	Pretty happy
9	Married	Good	117564.8	11.67	Very happy
11	Never married	Good	43560.0	10.68	Pretty happy

Interpretasi output: Delapan baris pertama data setelah transformasi. Kategori referensi yang dipilih adalah: Never married untuk marital dan Poor untuk health. Kolom log_income adalah transformasi logaritmik natural dari pendapatan riil; distribusinya jauh lebih simetris dibandingkan skala asli. Semua perbandingan odds ratio akan mengacu pada kategori-kategori referensi ini.

3.3.4 Distribusi Variabel Respons

data_bersih |>
  count(happy, name = "Jumlah") |>
  mutate(Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1)) |>
  rename(`Tingkat kebahagiaan` = happy) |>
  kable(caption = "Distribusi tingkat kebahagiaan responden") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE)

Distribusi tingkat kebahagiaan responden

Tingkat kebahagiaan	Jumlah	Proporsi
Not too happy	403	19.0%
Pretty happy	1242	58.5%
Very happy	479	22.6%

Interpretasi output: Distribusi tingkat kebahagiaan tidak merata: Pretty happy merupakan kategori terbanyak, diikuti Not too happy dan Very happy. Ketidakseimbangan ini lazim dalam data survei kebahagiaan dan tidak menghalangi penggunaan model ordinal.

3.4 Eksplorasi Data

3.4.1 Distribusi Tingkat Kebahagiaan

warna_happy <- c(
  "Not too happy" = "#5568b8",
  "Pretty happy"  = "#d18b2f",
  "Very happy"    = "#2f7f73"
)

data_bersih |>
  count(happy) |>
  mutate(proporsi = n / sum(n)) |>
  ggplot(aes(x = happy, y = proporsi, fill = happy)) +
  geom_col(width = 0.65, alpha = 0.94) +
  geom_text(aes(label = percent(proporsi, accuracy = 0.1)),
            vjust = -0.4, fontface = "bold", color = "#243142") +
  scale_y_continuous(labels = percent_format(), limits = c(0, 0.70)) +
  scale_fill_manual(values = warna_happy) +
  labs(title    = "Distribusi Tingkat Kebahagiaan",
       subtitle = "Respons ordinal: kategori memiliki urutan alami.",
       x = NULL, y = "Proporsi") +
  theme_minimal() +
  theme(legend.position = "none")

Interpretasi output: Mayoritas responden merasa cukup bahagia (Pretty happy). Kategori Very happy dan Not too happy masing-masing mewakili sekitar sepertiga dari sisa responden. Pola ini menunjukkan distribusi yang cenderung terpusat di kategori tengah.

3.4.2 Distribusi Kebahagiaan Menurut Status Pernikahan

prop_marital <- data_bersih |>
  count(marital, happy) |>
  group_by(marital) |>
  mutate(prop  = n / sum(n),
         label = percent(prop, accuracy = 1)) |>
  ungroup()

ggplot(prop_marital, aes(x = marital, y = prop, fill = happy)) +
  geom_col(width = 0.65, color = "white", linewidth = 0.6) +
  geom_text(aes(label = label),
            position = position_stack(vjust = 0.5),
            color = "white", fontface = "bold", size = 3.5) +
  scale_fill_manual(values = warna_happy) +
  scale_y_continuous(labels = percent_format()) +
  labs(title    = "Proporsi Tingkat Kebahagiaan Menurut Status Pernikahan",
       subtitle = "Responden menikah cenderung lebih bahagia.",
       x = "Status pernikahan", y = "Proporsi",
       fill = "Tingkat kebahagiaan") +
  theme_minimal() +
  theme(legend.position = "bottom")

Interpretasi output: Grafik stacked bar menampilkan perbedaan komposisi kebahagiaan antara responden menikah dan belum menikah. Responden yang menikah umumnya menunjukkan proporsi Very happy lebih tinggi dan Not too happy lebih rendah dibandingkan yang belum menikah, mengindikasikan bahwa status pernikahan adalah prediktor yang relevan.

3.4.3 Distribusi Kebahagiaan Menurut Kondisi Kesehatan

prop_health <- data_bersih |>
  count(health, happy) |>
  group_by(health) |>
  mutate(prop  = n / sum(n),
         label = ifelse(prop > 0.05, percent(prop, accuracy = 1), "")) |>
  ungroup()

ggplot(prop_health, aes(x = health, y = prop, fill = happy)) +
  geom_col(width = 0.70, color = "white", linewidth = 0.6) +
  geom_text(aes(label = label),
            position = position_stack(vjust = 0.5),
            color = "white", fontface = "bold", size = 3.2) +
  scale_fill_manual(values = warna_happy) +
  scale_y_continuous(labels = percent_format()) +
  labs(title    = "Proporsi Tingkat Kebahagiaan Menurut Kondisi Kesehatan",
       subtitle = "Kondisi kesehatan yang lebih baik berkaitan dengan kebahagiaan lebih tinggi.",
       x = "Kondisi kesehatan", y = "Proporsi",
       fill = "Tingkat kebahagiaan") +
  theme_minimal() +
  theme(legend.position = "bottom")

Interpretasi output: Terdapat gradien yang jelas: semakin baik kondisi kesehatan, semakin tinggi proporsi Very happy dan semakin rendah proporsi Not too happy. Pola ini konsisten dengan teori bahwa kesehatan merupakan determinan utama kebahagiaan subjektif.

3.4.4 Distribusi Log-Pendapatan Menurut Kebahagiaan

ggplot(data_bersih, aes(x = happy, y = log_income, fill = happy)) +
  geom_boxplot(alpha = 0.75, color = "#243142",
               outlier.alpha = 0.4, outlier.size = 1.2) +
  scale_fill_manual(values = warna_happy) +
  labs(title    = "Distribusi Log-Pendapatan Menurut Tingkat Kebahagiaan",
       subtitle = "Pendapatan lebih tinggi berkaitan dengan kebahagiaan lebih tinggi.",
       x = NULL, y = "Log pendapatan riil (ln USD)") +
  theme_minimal() +
  theme(legend.position = "none")

Interpretasi output: Boxplot menunjukkan bahwa median log-pendapatan cenderung meningkat seiring dengan tingginya kategori kebahagiaan. Responden Very happy memiliki median pendapatan lebih tinggi dibanding Not too happy. Hal ini mengindikasikan bahwa pendapatan memiliki hubungan positif dengan kebahagiaan, meskipun tidak selalu kuat.

3.4.5 Ringkasan Numerik Log-Pendapatan

data_bersih |>
  group_by(happy) |>
  summarise(
    Jumlah   = n(),
    Rata2    = mean(log_income, na.rm = TRUE),
    Median   = median(log_income, na.rm = TRUE),
    Std      = sd(log_income, na.rm = TRUE),
    Min      = min(log_income),
    Maks     = max(log_income),
    .groups  = "drop"
  ) |>
  rename(`Tingkat kebahagiaan` = happy) |>
  mutate(across(where(is.numeric), round, 3)) |>
  kable(caption = "Ringkasan log-pendapatan menurut tingkat kebahagiaan") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE)

Ringkasan log-pendapatan menurut tingkat kebahagiaan

Tingkat kebahagiaan	Jumlah	Rata2	Median	Std	Min	Maks
Not too happy	403	9.467	9.701	1.451	5.207	11.675
Pretty happy	1242	10.041	10.307	1.237	5.207	11.675
Very happy	479	10.137	10.307	1.307	5.207	11.675

Interpretasi output: Tabel memperlihatkan statistik deskriptif log-pendapatan untuk setiap kategori kebahagiaan. Perbedaan rata-rata dan median log-pendapatan antar kategori memberikan gambaran awal mengenai peran pendapatan sebagai prediktor.

3.5 Pembagian Training dan Testing

Data dibagi secara stratified 80:20 berdasarkan kategori happy agar proporsi setiap kategori terjaga pada kedua subset.

stratified_split <- function(y, prop = 0.8) {
  idx_by_class <- split(seq_along(y), y)
  train_idx <- lapply(idx_by_class,
    function(idx) sample(idx, size = floor(length(idx) * prop)))
  unlist(train_idx, use.names = FALSE)
}

set.seed(42)
train_id   <- stratified_split(as.integer(data_bersih$happy), prop = 0.8)
train_data <- data_bersih[ train_id, ]
test_data  <- data_bersih[-train_id, ]

bind_rows(
  train_data |> count(happy) |> mutate(Data = "Training"),
  test_data  |> count(happy) |> mutate(Data = "Testing")
) |>
  group_by(Data) |>
  mutate(Proporsi = scales::percent(n / sum(n), accuracy = 0.1)) |>
  ungroup() |>
  rename(`Tingkat kebahagiaan` = happy, Jumlah = n) |>
  kable(caption = "Distribusi kelas pada training dan testing") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE)

Distribusi kelas pada training dan testing

Tingkat kebahagiaan	Jumlah	Data	Proporsi
Not too happy	322	Training	19.0%
Pretty happy	993	Training	58.5%
Very happy	383	Training	22.6%
Not too happy	81	Testing	19.0%
Pretty happy	249	Testing	58.5%
Very happy	96	Testing	22.5%

Interpretasi output: Data training berisi 1698 observasi dan data testing berisi 426 observasi. Split stratified memastikan proporsi Not too happy, Pretty happy, dan Very happy pada training dan testing konsisten dengan proporsi data asli.

3.6 Model Regresi Logistik Ordinal

3.6.1 Persamaan Model untuk Studi Kasus Ini

Misalkan \(Y_i\) adalah tingkat kebahagiaan responden ke-\(i\), dengan:

\[ 1 = \text{Not too happy} < 2 = \text{Pretty happy} < 3 = \text{Very happy}. \]

Model cumulative logit membentuk dua persamaan (untuk \(J = 3\), ada \(J-1 = 2\) cutpoint):

\[ \operatorname{logit}\{P(Y_i \leq 1)\} = \alpha_1 + \beta_1 \cdot \text{log\_income}_i + \beta_2 \cdot \mathbf{1}[\text{marital}_i = \text{Married}] + \beta_3 \cdot \mathbf{1}[\text{health}_i = \text{Fair}] + \beta_4 \cdot \mathbf{1}[\text{health}_i = \text{Good}] + \beta_5 \cdot \mathbf{1}[\text{health}_i = \text{Excellent}], \]

\[ \operatorname{logit}\{P(Y_i \leq 2)\} = \alpha_2 + \beta_1 \cdot \text{log\_income}_i + \cdots + \beta_5 \cdot \mathbf{1}[\text{health}_i = \text{Excellent}]. \]

Karena proportional odds, koefisien \(\boldsymbol{\beta}\) sama untuk kedua persamaan; yang berbeda hanya cutpoint \(\alpha_1 < \alpha_2\).

3.6.2 Estimasi Model pada Data Training

fit_ord <- MASS::polr(
  happy ~ log_income + marital + health,
  data   = train_data,
  method = "logistic",
  Hess   = TRUE
)

# ── Ringkasan kecocokan model ──────────────────────────────────────────────
null_model <- MASS::polr(happy ~ 1,
                         data = train_data, method = "logistic", Hess = TRUE)

ll_null <- logLik(null_model)
ll_full <- logLik(fit_ord)
G2      <- -2 * (as.numeric(ll_null) - as.numeric(ll_full))
df_G2   <- attr(ll_full, "df") - attr(ll_null, "df")
p_G2    <- pchisq(G2, df = df_G2, lower.tail = FALSE)
pseudo_r2 <- 1 - as.numeric(ll_full) / as.numeric(ll_null)

data.frame(
  Keterangan = c("Jumlah observasi training",
                 "Log-likelihood model kosong",
                 "Log-likelihood model penuh",
                 "Statistik G² (uji rasio likelihood)",
                 "Derajat bebas G²",
                 "p-value uji G²",
                 "McFadden Pseudo-R²",
                 "AIC"),
  Nilai = c(
    nrow(train_data),
    round(as.numeric(ll_null), 3),
    round(as.numeric(ll_full), 3),
    round(G2, 3),
    df_G2,
    signif(p_G2, 3),
    round(pseudo_r2, 4),
    round(AIC(fit_ord), 3)
  )
) |>
  kable(caption = "Ringkasan kecocokan model regresi logistik ordinal") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE)

Ringkasan kecocokan model regresi logistik ordinal

Keterangan	Nilai
Jumlah observasi training	1698.0000
Log-likelihood model kosong	-1638.4480
Log-likelihood model penuh	-1522.2480
Statistik G² (uji rasio likelihood)	232.4000
Derajat bebas G²	5.0000
p-value uji G²	0.0000
McFadden Pseudo-R²	0.0709
AIC	3058.4950

Interpretasi output: Statistik \(G^2\) membandingkan model penuh dengan model kosong. Jika p-value < 0,05, prediktor secara bersama-sama berkontribusi signifikan dalam menjelaskan variasi tingkat kebahagiaan. McFadden Pseudo-\(R^2\) mengukur proporsi log-likelihood yang dijelaskan oleh prediktor; nilai di atas 0,10 sudah cukup substansial untuk data survei sosial.

3.6.3 Koefisien, Cutpoint, dan Odds Ratio

Model menghasilkan koefisien yang sama untuk kedua persamaan (proportional odds), ditambah dua cutpoint.

ord_coef <- coef(summary(fit_ord))

result_ord <- as.data.frame(ord_coef) |>
  tibble::rownames_to_column("parameter") |>
  dplyr::rename(estimate  = Value,
                std_error = `Std. Error`,
                z_value   = `t value`) |>
  mutate(
    p_value        = 2 * (1 - pnorm(abs(z_value))),
    jenis          = ifelse(grepl("\\|", parameter), "Cutpoint", "Koefisien"),
    estimate_slide = ifelse(jenis == "Koefisien", -estimate, estimate),
    OR_polr        = ifelse(jenis == "Koefisien", exp(estimate),   NA_real_),
    OR_slide       = ifelse(jenis == "Koefisien", exp(estimate_slide), NA_real_),
    CI_low_polr    = ifelse(jenis == "Koefisien",
                            exp(estimate - 1.96 * std_error), NA_real_),
    CI_high_polr   = ifelse(jenis == "Koefisien",
                            exp(estimate + 1.96 * std_error), NA_real_),
    Signifikan     = case_when(
      p_value < 0.001 ~ "***",
      p_value < 0.01  ~ "**",
      p_value < 0.05  ~ "*",
      p_value < 0.1   ~ ".",
      TRUE            ~ ""
    )
  )

result_ord |>
  mutate(across(c(estimate, estimate_slide, std_error, z_value,
                  p_value, OR_polr, OR_slide,
                  CI_low_polr, CI_high_polr), ~round(.x, 4))) |>
  kable(caption = "Ringkasan hasil regresi logistik ordinal (data training)",
        col.names = c("Parameter","Estimate polr","SE","z/t-value","p-value",
                      "Jenis","Estimate slide","OR polr","OR slide",
                      "CI polr 2.5%","CI polr 97.5%","Sig.")) |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = TRUE)

Ringkasan hasil regresi logistik ordinal (data training)

Parameter	Estimate polr	SE	z/t-value	p-value	Jenis	Estimate slide	OR polr	OR slide	CI polr 2.5%	CI polr 97.5%	Sig.
log_income	0.0322	0.0407	0.7933	0.4276	Koefisien	-0.0322	1.0328	0.9683	0.9537	1.1184
maritalMarried	1.0250	0.1091	9.3944	0.0000	Koefisien	-1.0250	2.7871	0.3588	2.2505	3.4516	***
healthFair	0.7007	0.3018	2.3220	0.0202	Koefisien	-0.7007	2.0152	0.4962	1.1154	3.6407
healthGood	1.4721	0.2965	4.9651	0.0000	Koefisien	-1.4721	4.3583	0.2294	2.4375	7.7927	***
healthExcellent	2.0995	0.3132	6.7028	0.0000	Koefisien	-2.0995	8.1618	0.1225	4.4174	15.0801	***
Not too happy|Pretty happy	0.6325	0.4573	1.3831	0.1666	Cutpoint	0.6325	NA	NA	NA	NA
Pretty happy|Very happy	3.6428	0.4671	7.7989	0.0000	Cutpoint	3.6428	NA	NA	NA	NA	***

Catatan: *** p < 0,001 | ** p < 0,01 | * p < 0,05 | . p < 0,1

Interpretasi kolom:

• Estimate polr: koefisien dari MASS::polr() dengan konvensi \(\alpha_j - \mathbf{x}^\top\boldsymbol{\beta}\).
• Estimate slide: tanda dibalik agar sesuai konvensi PDF/Agresti.
• OR polr: \(\exp(\hat{\beta}_{\text{polr}})\) — odds ratio dalam konvensi polr; OR > 1 berarti cenderung ke kategori lebih tinggi.
• OR slide: \(\exp(\hat{\beta}_{\text{slide}})\) — odds ratio dalam konvensi PDF; OR > 1 berarti cenderung ke kategori lebih rendah.
• CI polr: interval kepercayaan 95% dalam konvensi polr.

3.6.4 Cutpoint Model

result_ord |>
  filter(jenis == "Cutpoint") |>
  mutate(
    `Batas kumulatif` = gsub("\\|", " ≤ Y ≤ ", parameter),
    `Cutpoint (α)`    = round(estimate, 4),
    `SE`              = round(std_error, 4),
    `z-value`         = round(z_value, 4),
    `p-value`         = round(p_value, 4)
  ) |>
  select(`Batas kumulatif`, `Cutpoint (α)`, SE, `z-value`, `p-value`) |>
  kable(caption = "Cutpoint model ordinal") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE)

Cutpoint model ordinal

Batas kumulatif	Cutpoint (α)	SE	z-value	p-value
Not too happy ≤ Y ≤ Pretty happy	0.6325	0.4573	1.3831	0.1666
Pretty happy ≤ Y ≤ Very happy	3.6428	0.4671	7.7989	0.0000

Interpretasi output: Terdapat \(J - 1 = 2\) cutpoint karena respons memiliki 3 kategori. Cutpoint \(\alpha_1\) adalah batas antara Not too happy dan Pretty happy; \(\alpha_2\) adalah batas antara Pretty happy dan Very happy. Nilai \(\alpha_1 < \alpha_2\) harus terpenuhi agar model konsisten secara logis.

3.7 Prediksi dan Evaluasi

3.7.1 Prediksi Peluang pada Data Testing

prob_test  <- predict(fit_ord, newdata = test_data, type = "probs")
pred_kelas <- predict(fit_ord, newdata = test_data, type = "class")

head(
  data.frame(
    `Status aktual`       = test_data$happy,
    `Pred. kelas`         = pred_kelas,
    `P(Not too happy)`    = round(prob_test[,"Not too happy"], 4),
    `P(Pretty happy)`     = round(prob_test[,"Pretty happy"],  4),
    `P(Very happy)`       = round(prob_test[,"Very happy"],    4),
    check.names = FALSE
  ),
  8
) |>
  kable(caption = "Contoh prediksi peluang pada data testing") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE)

Contoh prediksi peluang pada data testing

Status aktual	Pred. kelas	P(Not too happy)	P(Pretty happy)	P(Very happy)
Pretty happy	Pretty happy	0.0961	0.5872	0.3167
Very happy	Pretty happy	0.0961	0.5872	0.3167
Pretty happy	Pretty happy	0.0537	0.4817	0.4646
Pretty happy	Pretty happy	0.0961	0.5872	0.3167
Pretty happy	Not too happy	0.5901	0.3768	0.0331
Not too happy	Pretty happy	0.4151	0.5199	0.0649
Very happy	Pretty happy	0.1870	0.6366	0.1765
Very happy	Pretty happy	0.3997	0.5314	0.0689

Interpretasi output: Setiap baris menampilkan tiga peluang prediksi yang jumlahnya selalu = 1. Model mengklasifikasikan observasi ke kategori dengan peluang tertinggi. Pada respons ordinal, peluang yang tersebar di antara kategori berdekatan lebih informatif daripada sekadar kelas prediksi.

3.7.2 Confusion Matrix

cm <- addmargins(table(
  `Aktual`   = test_data$happy,
  `Prediksi` = pred_kelas
))

cm |>
  kable(caption = "Confusion matrix data testing (baris = aktual, kolom = prediksi)") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed","responsive"),
                full_width = FALSE)

Confusion matrix data testing (baris = aktual, kolom = prediksi)

	Not too happy	Pretty happy	Very happy	Sum
Not too happy	2	79	0	81
Pretty happy	5	244	0	249
Very happy	1	95	0	96
Sum	8	418	0	426

Interpretasi output: Diagonal confusion matrix menunjukkan prediksi yang benar untuk setiap kategori. Pada respons ordinal, kesalahan prediksi tidak selalu sama beratnya: memprediksi Very happy padahal aktual Not too happy jauh lebih besar kesalahannya dibandingkan memprediksi Pretty happy padahal aktual Not too happy.

3.7.3 Metrik Evaluasi

pred_kelas <- predict(fit_ord, newdata = test_data, type = "class")
actual_vec <- as.character(test_data$happy)
pred_vec   <- as.character(pred_kelas)

overall_acc  <- mean(actual_vec == pred_vec)
kelas_levels <- c("Not too happy", "Pretty happy", "Very happy")

per_kelas <- lapply(kelas_levels, function(k) {
  tp <- sum(actual_vec == k & pred_vec == k)
  fp <- sum(actual_vec != k & pred_vec == k)
  fn <- sum(actual_vec == k & pred_vec != k)
  tn <- sum(actual_vec != k & pred_vec != k)
  sensitivity <- ifelse((tp + fn) == 0, NA_real_, tp / (tp + fn))
  precision   <- ifelse((tp + fp) == 0, NA_real_, tp / (tp + fp))
  specificity <- ifelse((tn + fp) == 0, NA_real_, tn / (tn + fp))
  f1          <- ifelse(is.na(precision) | is.na(sensitivity) |
                          (precision + sensitivity) == 0, NA_real_,
                        2 * precision * sensitivity /
                          (precision + sensitivity))
  data.frame(Kategori    = k,
             Sensitivity = round(sensitivity, 3),
             Specificity = round(specificity, 3),
             Precision   = round(precision, 3),
             F1          = round(f1, 3))
}) |> bind_rows()

knitr::kable(per_kelas,
             caption = "Metrik evaluasi per kategori (data testing)")

Metrik evaluasi per kategori (data testing)

Kategori	Sensitivity	Specificity	Precision	F1
Not too happy	0.025	0.983	0.250	0.045
Pretty happy	0.980	0.017	0.584	0.732
Very happy	0.000	1.000	NA	NA

score_map <- setNames(1:3, kelas_levels)
mae_ord   <- mean(abs(score_map[actual_vec] - score_map[pred_vec]),
                  na.rm = TRUE)

knitr::kable(
  data.frame(
    Metrik = c("Overall Accuracy",
               "Macro-avg Sensitivity",
               "Macro-avg Precision",
               "Macro-avg F1-score",
               "Mean Absolute Error (ordinal)"),
    Nilai  = round(c(
      overall_acc,
      mean(per_kelas$Sensitivity, na.rm = TRUE),
      mean(per_kelas$Precision,   na.rm = TRUE),
      mean(per_kelas$F1,          na.rm = TRUE),
      mae_ord
    ), 3)
  ),
  caption = "Metrik evaluasi keseluruhan"
)

Metrik evaluasi keseluruhan

Metrik	Nilai
Overall Accuracy	0.577
Macro-avg Sensitivity	0.335
Macro-avg Precision	0.417
Macro-avg F1-score	0.388
Mean Absolute Error (ordinal)	0.425

Interpretasi output — metrik per kategori: Sensitivity mengukur kemampuan model mengenali observasi yang benar-benar masuk kategori tersebut. Precision mengukur ketepatan prediksi positif. F1-score menyeimbangkan keduanya.

Interpretasi output — MAE ordinal: Mean Absolute Error pada skor ordinal mengukur rata-rata selisih antara posisi kategori aktual dan prediksi. MAE = 0 berarti sempurna; MAE = 1 berarti rata-rata salah satu tingkat kategori. Metrik ini lebih sesuai untuk respons ordinal dibandingkan akurasi sederhana.

3.7.4 Visualisasi Confusion Matrix

as.data.frame(table(
  Aktual   = test_data$happy,
  Prediksi = pred_kelas
)) |>
  group_by(Aktual) |>
  mutate(prop = Freq / sum(Freq)) |>
  ungroup() |>
  ggplot(aes(x = Prediksi, y = Aktual, fill = prop)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(aes(label = paste0(Freq, "\n(", percent(prop, accuracy = 1), ")")),
            fontface = "bold", size = 3.4) +
  scale_fill_gradient(low = "#f8f9fa", high = "#2f7f73",
                      labels = percent_format()) +
  labs(title    = "Confusion Matrix (Proporsi per Baris Aktual)",
       subtitle = "Angka = frekuensi; persentase = dari total aktual per kategori.",
       x = "Prediksi", y = "Aktual", fill = "Proporsi") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))

Interpretasi output: Sel diagonal yang lebih gelap menandakan prediksi yang lebih akurat. Perhatikan apakah kesalahan prediksi cenderung terjadi antara kategori berdekatan (misalnya Pretty happy diprediksi Very happy) atau antara kategori jauh (Not too happy diprediksi Very happy). Pada model ordinal yang baik, kesalahan seharusnya lebih banyak terjadi pada kategori berdekatan.

3.8 Visualisasi Probabilitas Prediksi

3.8.1 Peluang Kategori terhadap Log-Pendapatan

grid_income <- expand.grid(
  log_income = seq(min(data_bersih$log_income),
                   max(data_bersih$log_income),
                   length.out = 120),
  marital    = "Married",
  health     = "Good"
)

grid_prob_income <- predict(fit_ord, newdata = grid_income, type = "probs")

grid_income |>
  bind_cols(as.data.frame(grid_prob_income)) |>
  tidyr::pivot_longer(
    cols      = c("Not too happy", "Pretty happy", "Very happy"),
    names_to  = "kategori",
    values_to = "probabilitas"
  ) |>
  mutate(kategori = factor(kategori,
                           levels = c("Not too happy","Pretty happy","Very happy"))) |>
  ggplot(aes(x = log_income, y = probabilitas, color = kategori)) +
  geom_line(linewidth = 1.25) +
  scale_y_continuous(labels = percent_format()) +
  scale_color_manual(values = warna_happy) +
  labs(title    = "Prediksi Probabilitas Tingkat Kebahagiaan vs Log-Pendapatan",
       subtitle = "Marital = Married; Health = Good. Variabel lain ditahan konstan.",
       x = "Log pendapatan riil (ln USD)",
       y = "Probabilitas prediksi",
       color = "Tingkat kebahagiaan") +
  theme_minimal()

Interpretasi output: Grafik ini menampilkan \(P(Y = j)\) sebagai fungsi log-pendapatan. Jika garis Very happy meningkat dan Not too happy menurun seiring kenaikan pendapatan, model menunjukkan bahwa pendapatan lebih tinggi dikaitkan dengan kebahagiaan lebih tinggi. Kurva kategori tunggal tidak harus sejajar — yang paralel adalah cumulative logit pada skala logit, bukan probabilitas kategori \(P(Y = j)\).

3.8.2 Parallel Lines: Cumulative Logit

eps <- 1e-6

grid_income |>
  bind_cols(as.data.frame(grid_prob_income)) |>
  rename(p1 = `Not too happy`, p2 = `Pretty happy`, p3 = `Very happy`) |>
  mutate(
    `P(Y ≤ Not too happy)` = p1,
    `P(Y ≤ Pretty happy)`  = p1 + p2
  ) |>
  tidyr::pivot_longer(
    cols      = c(`P(Y ≤ Not too happy)`, `P(Y ≤ Pretty happy)`),
    names_to  = "batas_kumulatif",
    values_to = "prob_kumulatif"
  ) |>
  mutate(
    prob_kumulatif  = pmin(pmax(prob_kumulatif, eps), 1 - eps),
    cumulative_logit = qlogis(prob_kumulatif)
  ) |>
  ggplot(aes(x = log_income, y = cumulative_logit,
             color = batas_kumulatif)) +
  geom_line(linewidth = 1.25) +
  scale_color_manual(values = c(
    "P(Y ≤ Not too happy)" = "#2f7f73",
    "P(Y ≤ Pretty happy)"  = "#d18b2f"
  )) +
  labs(title    = "Parallel Lines pada Model Ordinal",
       subtitle = "Garis yang sejajar adalah cumulative logit, bukan probabilitas kategori tunggal.",
       x = "Log pendapatan riil (ln USD)",
       y = "Logit peluang kumulatif",
       color = "Batas kumulatif") +
  theme_minimal()

Interpretasi output: Dua garis yang sejajar pada grafik ini merupakan bukti visual asumsi proportional odds terpenuhi. Kedua garis memiliki slope yang sama — hanya intercept (\(\alpha_j\)) yang berbeda. Jika dua garis tersebut tidak sejajar secara substansial, asumsi proportional odds patut dipertanyakan.

3.9 Pemeriksaan Asumsi Proportional Odds

3.9.1 Definisi Asumsi

Asumsi utama model ordinal adalah proportional odds atau parallel lines: efek prediktor dianggap sama untuk setiap batas kumulatif. Untuk tiga kategori kebahagiaan, model membentuk dua batas:

\[ Y \leq \text{Not too happy} \quad \text{dan} \quad Y \leq \text{Pretty happy}. \]

Asumsi proportional odds menyatakan bahwa slope \(\boldsymbol{\beta}\) sama untuk kedua batas tersebut.

3.9.2 Uji dengan Paket ordinal

fit_clm <- ordinal::clm(
  happy ~ log_income + marital + health,
  data = train_data,
  link = "logit"
)

# Uji indikatif apakah efek prediktor perlu dibuat non-proportional
nominal_test_result <- ordinal::nominal_test(fit_clm)
print(nominal_test_result)

## Tests of nominal effects
## 
## formula: happy ~ log_income + marital + health
##            Df  logLik    AIC     LRT  Pr(>Chi)    
## <none>        -1522.2 3058.5                      
## log_income  1 -1515.0 3045.9 14.5854 0.0001339 ***
## marital     1 -1522.1 3060.3  0.2385 0.6252956    
## health      3 -1515.3 3050.5 13.9645 0.0029539 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretasi output: Fungsi nominal_test() menguji apakah koefisien setiap prediktor perlu dibuat berbeda untuk setiap batas kumulatif (non-proportional). p-value kecil (< 0,05) untuk suatu prediktor menjadi indikasi bahwa efek prediktor tersebut mungkin tidak sepenuhnya memenuhi asumsi proportional odds.

Namun hasil uji tidak boleh dibaca secara mekanis. Keputusan perlu mempertimbangkan: (a) teori substantif, (b) ukuran sampel yang besar cenderung membuat uji sangat sensitif, (c) visualisasi cumulative logit, dan (d) tujuan analisis.

3.9.3 Jika Asumsi Terlanggar

Jika asumsi proportional odds tidak terpenuhi, alternatif yang dapat dipertimbangkan:

• Partial proportional odds model — efek prediktor tertentu dibiarkan berbeda antar batas.
• Generalized ordered logit — semua koefisien berbeda antar batas (setara model multinomial pada data ordinal).
• Adjacent-category logit — membandingkan kategori berdekatan.
• Continuation-ratio logit — cocok untuk proses bertahap/sekuensial.
• Multinomial logistic regression — sebagai pembanding non-ordinal (tidak mengeksploitasi urutan).

3.10 Interpretasi Hasil

3.10.1 Kebaikan Model

Berdasarkan hasil estimasi pada data training (n = 1698), model secara keseluruhan signifikan berdasarkan uji rasio likelihood (\(G^2 =\) 232.4; df = 5; p = 3.27^{-48}). Ketiga prediktor — log-pendapatan, status pernikahan, dan kondisi kesehatan — secara bersama-sama memberikan kontribusi nyata dalam menjelaskan variasi tingkat kebahagiaan.

McFadden Pseudo-\(R^2\) sebesar 0.071 menunjukkan proporsi log-likelihood yang berhasil dijelaskan model. Nilai ini wajar untuk data survei sosial; kebahagiaan subjektif dipengaruhi banyak faktor yang tidak tersedia dalam dataset.

3.10.2 Koefisien: Log-Pendapatan

Koefisien log_income pada output polr() positif berarti (dalam konvensi polr) kenaikan log-pendapatan meningkatkan odds kumulatif \(P(Y \leq j)\), sehingga responden cenderung bergeser ke kategori kebahagiaan lebih tinggi. Dalam konvensi slide, tandanya terbalik namun interpretasi arah tetap sama.

Setiap kenaikan satu satuan log-pendapatan mengalikan odds kumulatif sebesar \(\exp(\hat{\beta}_{\text{polr}})\). Jika OR > 1 (polr), maka pendapatan yang lebih tinggi dikaitkan dengan kecenderungan lebih kuat untuk berada pada kategori kebahagiaan yang lebih tinggi.

3.10.3 Koefisien: Status Pernikahan

Koefisien maritalMarried yang positif pada output polr menunjukkan bahwa responden yang menikah memiliki odds kumulatif \(P(Y \leq j)\) lebih tinggi dibandingkan yang tidak menikah (referensi). Artinya, dalam konvensi polr, menikah dikaitkan dengan kecenderungan lebih besar berada di kategori kebahagiaan lebih tinggi.

3.10.4 Koefisien: Kondisi Kesehatan

Level-level health menunjukkan gradien yang konsisten: semakin baik kondisi kesehatan (dari Poor sebagai referensi menuju Excellent), semakin besar kecenderungan berada di kategori kebahagiaan lebih tinggi. Level dengan p-value < 0,05 memberikan bukti statistik yang kuat. Kondisi kesehatan kemungkinan merupakan prediktor terkuat dalam model ini mengingat pola gradien yang jelas pada eksplorasi data.

3.10.5 Evaluasi Klasifikasi

Overall accuracy pada data testing adalah 57.7%. MAE ordinal sebesar 0.425 berarti rata-rata prediksi model meleset sekitar 0.425 tingkat kategori dari nilai aktual. Nilai MAE mendekati 0 menunjukkan prediksi yang rata-rata cukup dekat dengan kategori yang benar.

Performa model tidak merata antar kategori: Pretty happy (kategori terbanyak) umumnya diprediksi paling baik, sementara Not too happy (minoritas) cenderung lebih sulit dideteksi. Kesalahan prediksi diharapkan lebih banyak terjadi antara kategori berdekatan (misalnya Not too happy ↔︎ Pretty happy) daripada antara kategori berjauhan (Not too happy ↔︎ Very happy).

3.11 Kesimpulan

Analisis regresi logistik ordinal terhadap data GSS 2024 (n = 2124) menghasilkan beberapa temuan utama:

Pertama, signifikansi model. Model secara keseluruhan signifikan (\(G^2 =\) 232.4; df = 5; p = 3.27^{-48}), menunjukkan bahwa log-pendapatan, status pernikahan, dan kondisi kesehatan secara bersama-sama berpengaruh terhadap tingkat kebahagiaan.

Kedua, prediktor dominan. Kondisi kesehatan kemungkinan merupakan prediktor terkuat berdasarkan gradien yang konsisten pada eksplorasi data: semakin baik kesehatan, semakin tinggi kecenderungan Very happy. Status pernikahan juga signifikan: responden yang menikah cenderung lebih bahagia. Log-pendapatan memberikan kontribusi positif, konsisten dengan teori utilitas marginal yang semakin menurun.

Ketiga, performa klasifikasi. Model mencapai overall accuracy 57.7% pada data testing, dengan MAE ordinal 0.425. Ketimpangan performa antar kategori mencerminkan ketidakseimbangan distribusi kelas.

Keempat, asumsi dan keterbatasan. Asumsi proportional odds perlu diperiksa secara formal. McFadden Pseudo-\(R^2\) sebesar 0.071 menunjukkan masih banyak variasi yang belum dijelaskan, yang wajar mengingat kompleksitas kebahagiaan subjektif. Model ini memberikan gambaran yang informatif dan parsimonious mengenai determinan demografis kebahagiaan, namun tidak boleh diinterpretasikan sebagai hubungan kausal.

4 Regresi Logistik Poisson

Regresi Poisson adalah model regresi untuk variabel respons berupa data hitung (count data), yaitu respons yang bernilai bilangan bulat nonnegatif: \[Y_i \in \{0, 1, 2, \ldots\}.\]

Contoh respons yang cocok dianalisis dengan regresi Poisson: jumlah anggota rumah tangga, jumlah kunjungan pasien ke klinik, jumlah kecelakaan di suatu ruas jalan, jumlah klaim asuransi, atau jumlah komplain pelanggan.

Model Poisson digunakan ketika kita ingin menjelaskan rata-rata jumlah kejadian sebagai fungsi dari beberapa prediktor. Berbeda dari regresi logistik yang memodelkan peluang kategori, regresi Poisson memodelkan nilai harapan hitungan: \[E(Y_i \mid \mathbf{x}_i) = \mu_i.\]

Distribusi dan Model

Asumsi dasar model Poisson adalah: \[Y_i \mid \mathbf{x}_i \sim \text{Poisson}(\mu_i),\] dengan fungsi peluang: \[P(Y_i = y_i \mid \mathbf{x}_i) = \frac{\exp(-\mu_i)\,\mu_i^{y_i}}{y_i!}, \quad y_i = 0, 1, 2, \ldots.\]

Pada distribusi Poisson berlaku sifat equidispersion, yaitu rata-rata dan varians kondisional bernilai sama: \[E(Y_i \mid \mathbf{x}_i) = \mu_i, \qquad \text{Var}(Y_i \mid \mathbf{x}_i) = \mu_i.\]

Regresi Poisson menggunakan fungsi penghubung log: \[\log(\mu_i) = \eta_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}.\]

Dengan demikian: \[\mu_i = \exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}).\]

Fungsi log dipakai karena \(\mu_i\) harus selalu positif. Dengan bentuk eksponensial, prediksi rata-rata hitungan tidak mungkin negatif.

Model dengan Exposure dan Offset

Dalam banyak kasus, jumlah kejadian bergantung pada lama waktu pengamatan atau ukuran populasi yang berisiko. Misalnya: jumlah kecelakaan per tahun pengamatan, jumlah klaim per jumlah polis, atau jumlah kasus penyakit per jumlah penduduk.

Jika \(t_i\) adalah exposure, model Poisson ditulis: \[\mu_i = t_i \lambda_i,\] dengan \(\lambda_i\) adalah rate kejadian per satu unit exposure. Karena \(\log(\mu_i) = \log(t_i) + \log(\lambda_i)\), maka model menjadi: \[\log(\mu_i) = \log(t_i) + \beta_0 + \mathbf{x}_i^\top\boldsymbol{\beta}.\]

Komponen \(\log(t_i)\) disebut offset, yaitu variabel yang koefisiennya dipaksa sama dengan 1. Pada analisis ini tidak digunakan offset karena seluruh observasi berasal dari survei lintas-seksi yang diasumsikan memiliki unit eksposur yang setara.

Studi kasus: Analisis faktor ekonomi dan latar belakang budaya yang memengaruhi jumlah orang yang tinggal dalam satu rumah tangga berdasarkan data General Social Survey (GSS) 2024.

Tujuan analisis: Membangun model regresi Poisson untuk menjelaskan rata-rata jumlah anggota rumah tangga (hompop) sebagai fungsi dari pendapatan keluarga (realinc) dan ras (race).

4.1 Ringkasan Masalah

Ukuran rumah tangga merupakan variabel hitungan (count) yang bernilai bilangan bulat positif, sehingga cocok dianalisis dengan regresi Poisson. Variabel respons pada analisis ini adalah jumlah orang yang tinggal dalam satu rumah tangga (hompop): \[Y_i = \text{jumlah anggota rumah tangga ke-}i \in \{1, 2, 3, \ldots\}.\]

Model yang digunakan: \[\log(\mu_i) = \beta_0 + \beta_1 \log(\texttt{realinc}_i) + \beta_2 \texttt{raceWhite}_i,\] dengan \(\mu_i = E(\texttt{hompop}_i \mid \mathbf{x}_i)\) adalah rata-rata jumlah anggota rumah tangga untuk responden ke-\(i\).

4.2 Sumber Data

Dataset berasal dari General Social Survey (GSS) 2024, survei sosial nasional Amerika Serikat yang diselenggarakan oleh NORC di University of Chicago. Dataset terdiri dari 1.289 responden. Dari lima variabel yang tersedia, analisis ini menggunakan tiga variabel: hompop, race, dan realinc.

4.3 Persiapan Data

4.3.1 Membaca Data dari Excel

raw_data <- readxl::read_excel("Data Poisson.xlsx", sheet = 1)

ringkasan_data <- data.frame(
  Keterangan = c("Jumlah observasi", "Jumlah variabel yang digunakan"),
  Nilai      = c(nrow(raw_data), 3)
)

knitr::kable(ringkasan_data, caption = "Ukuran dataset GSS 2024 (Poisson)")

Ukuran dataset GSS 2024 (Poisson)

Keterangan	Nilai
Jumlah observasi	1289
Jumlah variabel yang digunakan	3

Interpretasi output: Dataset berisi 1289 observasi. Dari seluruh variabel yang tersedia, analisis ini hanya menggunakan hompop sebagai respons, race dan realinc sebagai prediktor.

4.3.2 Penjelasan Variabel

kamus <- data.frame(
  `Nama kolom`          = c("race", "hompop", "realinc"),
  `Keterangan`          = c(
    "Ras responden",
    "Jumlah orang dalam satu rumah tangga (variabel respons)",
    "Pendapatan keluarga dalam setahun (USD)"
  ),
  `Kategori / Rentang`  = c(
    "Black, White",
    "1 – 10 orang",
    "181,5 – 117.564,8 USD"
  ),
  `Tipe dalam analisis` = c(
    "Kategorik (2 level); Black sebagai referensi",
    "Respons hitung (count)",
    "Numerik kontinu; ditransformasi log"
  ),
  check.names = FALSE
)

knitr::kable(kamus, caption = "Penjelasan variabel yang digunakan dalam analisis")

Penjelasan variabel yang digunakan dalam analisis

Nama kolom	Keterangan	Kategori / Rentang	Tipe dalam analisis
race	Ras responden	Black, White	Kategorik (2 level); Black sebagai referensi
hompop	Jumlah orang dalam satu rumah tangga (variabel respons)	1 – 10 orang	Respons hitung (count)
realinc	Pendapatan keluarga dalam setahun (USD)	181,5 – 117.564,8 USD	Numerik kontinu; ditransformasi log

Interpretasi output: Terdapat satu variabel respons hitung (hompop) dan dua prediktor: satu numerik (realinc) dan satu kategorik (race). Variabel realinc ditransformasi logaritma untuk menstabilkan distribusinya yang sangat miring ke kanan (right-skewed).

4.3.3 Transformasi dan Pembersihan Data

data_bersih <- raw_data |>
  select(id_, hompop, race, realinc) |>
  mutate(
    hompop  = as.integer(hompop),
    realinc = as.numeric(realinc),
    log_inc = log(realinc),
    race    = factor(race, levels = c("Black", "White"))
  ) |>
  tidyr::drop_na()

knitr::kable(
  head(data_bersih |> select(id_, hompop, race, realinc, log_inc), 8),
  caption = "Contoh data setelah transformasi",
  digits  = 2
)

Contoh data setelah transformasi

id_	hompop	race	realinc	log_inc
1	3	Black	43560.0	10.68
3	2	White	43560.0	10.68
5	1	White	24502.5	10.11
8	2	White	3267.0	8.09
12	1	White	3267.0	8.09
16	3	White	36300.0	10.50
19	2	White	43560.0	10.68
20	4	White	24502.5	10.11

Interpretasi output: Delapan baris pertama data setelah transformasi. Kategori referensi yang dipilih adalah Black untuk race. Kolom log_inc merupakan logaritma natural pendapatan yang digunakan sebagai prediktor numerik dalam model.

4.4 Distribusi Respons Hitungan

dist_hompop <- data_bersih |>
  count(hompop, name = "Frekuensi") |>
  mutate(Proporsi = scales::percent(Frekuensi / sum(Frekuensi),
                                    accuracy = 0.1)) |>
  rename(`Jumlah anggota RT` = hompop)

knitr::kable(dist_hompop,
             caption = "Distribusi jumlah anggota rumah tangga")

Distribusi jumlah anggota rumah tangga

Jumlah anggota RT	Frekuensi	Proporsi
1	457	35.5%
2	442	34.3%
3	162	12.6%
4	136	10.6%
5	49	3.8%
6	25	1.9%
7	12	0.9%
8	4	0.3%
9	1	0.1%
10	1	0.1%

Interpretasi output: Distribusi hompop menunjukkan pola khas data hitung: nilai kecil mendominasi (1–3 orang) dan frekuensi menurun tajam untuk nilai yang lebih besar. Pola ini konsisten dengan asumsi distribusi Poisson.

ggplot(data_bersih, aes(x = factor(hompop))) +
  geom_bar(fill = "#2a9d8f", color = "white", linewidth = 0.7,
           alpha = 0.92) +
  geom_text(stat = "count",
            aes(label = after_stat(count)),
            vjust = -0.4, fontface = "bold", color = "#243142") +
  labs(title    = "Distribusi Jumlah Anggota Rumah Tangga",
       subtitle = "Respons Poisson berupa hitungan bilangan bulat positif.",
       x = "Jumlah anggota rumah tangga",
       y = "Frekuensi") +
  theme_minimal(base_size = 12)

Interpretasi output: Grafik memperlihatkan distribusi data hitung yang condong ke kanan, dengan nilai 1 dan 2 paling sering muncul. Tidak terdapat nilai nol karena setiap responden setidaknya tinggal sendirian (hompop ≥ 1).

4.5 Eksplorasi Data

4.5.1 Jumlah Anggota RT Menurut Ras

ggplot(data_bersih, aes(x = race, y = hompop, fill = race)) +
  geom_boxplot(alpha = 0.75, color = "#243142",
               outlier.alpha = 0.4, outlier.size = 1.2) +
  scale_fill_manual(values = c("Black" = "#e76f51", "White" = "#2a9d8f")) +
  labs(title = "Distribusi Jumlah Anggota RT Menurut Ras",
       x = NULL, y = "Jumlah anggota rumah tangga") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

Interpretasi output: Boxplot menampilkan perbedaan sebaran hompop antara kelompok Black dan White. Perbedaan median atau rentang interkuartil antar kelompok memberikan gambaran awal apakah ras berkontribusi pada ukuran rumah tangga.

4.5.2 Hubungan Pendapatan dan Jumlah Anggota RT

ggplot(data_bersih, aes(x = log_inc, y = hompop, color = race)) +
  geom_point(alpha = 0.35, size = 1.4) +
  geom_smooth(method = "loess", se = FALSE, linewidth = 1.2) +
  scale_color_manual(values = c("Black" = "#e76f51", "White" = "#2a9d8f")) +
  labs(title    = "Hubungan Log-Pendapatan dan Jumlah Anggota RT",
       subtitle = "Garis LOESS menunjukkan pola hubungan per kelompok ras.",
       x = "Log pendapatan (ln USD)",
       y = "Jumlah anggota rumah tangga",
       color = "Ras") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Interpretasi output: Scatter plot dengan kurva LOESS menampilkan pola umum hubungan antara log-pendapatan dan jumlah anggota rumah tangga untuk masing-masing kelompok ras. Jika kurva cenderung turun, maka pendapatan yang lebih tinggi dikaitkan dengan rumah tangga yang lebih kecil.

4.5.3 Ringkasan Numerik Menurut Ras

ringkasan_ras <- data_bersih |>
  group_by(race) |>
  summarise(
    Jumlah  = n(),
    Rata2   = mean(hompop),
    Median  = median(hompop),
    Varians = var(hompop),
    Std     = sd(hompop),
    Min     = min(hompop),
    Maks    = max(hompop),
    .groups = "drop"
  ) |>
  mutate(across(where(is.numeric), round, 2))

knitr::kable(ringkasan_ras,
             caption = "Ringkasan jumlah anggota RT menurut ras")

Ringkasan jumlah anggota RT menurut ras

race	Jumlah	Rata2	Median	Varians	Std	Min	Maks
Black	264	2.06	2	1.87	1.37	1	7
White	1025	2.30	2	1.93	1.39	1	10

Interpretasi output: Tabel merangkum statistik deskriptif hompop per kelompok ras. Pada model Poisson, asumsi equidispersion menyatakan rata-rata dan varians kondisional seharusnya serupa. Jika varians jauh lebih besar dari rata-rata, ini menjadi indikasi awal overdispersion.

4.5.4 Distribusi Pendapatan Menurut Ras

ggplot(data_bersih, aes(x = race, y = log_inc, fill = race)) +
  geom_violin(alpha = 0.6, color = "#243142", trim = FALSE) +
  geom_boxplot(width = 0.12, fill = "white",
               color = "#243142", outlier.size = 1) +
  scale_fill_manual(values = c("Black" = "#e76f51", "White" = "#2a9d8f")) +
  labs(title = "Distribusi Log-Pendapatan Menurut Ras",
       x = NULL, y = "Log pendapatan (ln USD)") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

Interpretasi output: Violin plot menampilkan sebaran log-pendapatan untuk masing-masing kelompok ras. Perbedaan distribusi pendapatan antar ras perlu dicermati karena dapat berinteraksi dengan ukuran rumah tangga dalam model.

4.6 Estimasi Model Poisson dengan R

Model diestimasi menggunakan fungsi glm() dengan family = poisson. Pendapatan digunakan dalam skala log. Tidak ada offset karena seluruh observasi berasal dari survei lintas-seksi yang diasumsikan memiliki eksposur setara.

model_poisson <- glm(
  hompop ~ log_inc + race,
  data   = data_bersih,
  family = poisson(link = "log")
)

summary(model_poisson)

## 
## Call:
## glm(formula = hompop ~ log_inc + race, family = poisson(link = "log"), 
##     data = data_bersih)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.07960    0.16219  -0.491    0.624    
## log_inc      0.08602    0.01664   5.170 2.33e-07 ***
## raceWhite    0.04272    0.04921   0.868    0.385    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 966.86  on 1288  degrees of freedom
## Residual deviance: 933.72  on 1286  degrees of freedom
## AIC: 4269.4
## 
## Number of Fisher Scoring iterations: 5

Interpretasi output: glm() menampilkan koefisien, standard error, z-value, dan p-value untuk setiap prediktor. Baris Residual deviance dan Null deviance digunakan untuk menilai kontribusi model secara keseluruhan.

4.7 Ringkasan Koefisien dan IRR

coef_pois <- as.data.frame(coef(summary(model_poisson))) |>
  tibble::rownames_to_column("Parameter") |>
  dplyr::rename(
    Estimate = Estimate,
    SE       = `Std. Error`,
    z_value  = `z value`,
    p_value  = `Pr(>|z|)`
  ) |>
  mutate(
    IRR           = exp(Estimate),
    CI_low        = exp(Estimate - 1.96 * SE),
    CI_high       = exp(Estimate + 1.96 * SE),
    Perubahan_pct = 100 * (IRR - 1),
    Signifikan    = case_when(
      p_value < 0.001 ~ "***",
      p_value < 0.01  ~ "**",
      p_value < 0.05  ~ "*",
      p_value < 0.1   ~ ".",
      TRUE            ~ ""
    )
  )

knitr::kable(
  coef_pois |>
    mutate(across(c(Estimate, SE, z_value, p_value,
                    IRR, CI_low, CI_high, Perubahan_pct),
                  round, 4)),
  caption   = "Ringkasan hasil regresi Poisson",
  col.names = c("Parameter", "Estimate", "SE", "z-value", "p-value",
                "IRR", "CI 2,5%", "CI 97,5%", "Perubahan (%)", "Sig.")
)

Ringkasan hasil regresi Poisson

Parameter	Estimate	SE	z-value	p-value	IRR	CI 2,5%	CI 97,5%	Perubahan (%)	Sig.
(Intercept)	-0.0796	0.1622	-0.4908	0.6236	0.9235	0.6720	1.2691	-7.6513
log_inc	0.0860	0.0166	5.1705	0.0000	1.0898	1.0549	1.1259	8.9827	***
raceWhite	0.0427	0.0492	0.8681	0.3853	1.0437	0.9477	1.1493	4.3650

Catatan: *** p < 0,001 | ** p < 0,01 | * p < 0,05 | . p < 0,1

Interpretasi output: Kolom IRR adalah nilai eksponensial koefisien. IRR > 1 berarti prediktor meningkatkan rata-rata hompop; IRR < 1 berarti prediktor menurunkannya. Kolom Perubahan (%) menyatakan perubahan persentase pada rata-rata hitungan untuk setiap kenaikan satu unit prediktor atau perpindahan ke kategori non-referensi.

4.8 Interpretasi Koefisien

Log-pendapatan (log_inc): Setiap kenaikan satu unit log-pendapatan dikaitkan dengan perubahan rata-rata hompop sebesar \(\{(\text{IRR} - 1) \times 100\}\%\), dengan ras konstan. Jika koefisiennya negatif, maka rumah tangga berpendapatan lebih tinggi cenderung memiliki lebih sedikit anggota — kemungkinan mencerminkan perbedaan gaya hidup atau biaya hidup yang lebih tinggi di wilayah perkotaan.

Ras (raceWhite): Koefisien ini membandingkan rata-rata hompop kelompok White terhadap Black sebagai kategori referensi, setelah mengendalikan pendapatan. IRR < 1 berarti rata-rata ukuran rumah tangga White lebih kecil daripada Black; IRR > 1 berarti sebaliknya.

4.8.1 Visualisasi IRR

irr_plot <- coef_pois |>
  filter(Parameter != "(Intercept)")

ggplot(irr_plot,
       aes(x = IRR, y = reorder(Parameter, IRR))) +
  geom_vline(xintercept = 1, linetype = "dashed",
             color = "#6c757d", linewidth = 0.8) +
  geom_errorbarh(aes(xmin = CI_low, xmax = CI_high),
                 height = 0.2, linewidth = 0.8,
                 color = "#2a9d8f") +
  geom_point(size = 4, color = "#2a9d8f") +
  scale_x_log10() +
  labs(title    = "Incidence Rate Ratio (IRR) Model Poisson",
       subtitle = "Skala logaritmik; garis putus-putus = IRR 1 (tidak ada efek).",
       x        = "IRR (skala log)",
       y        = NULL) +
  theme_minimal(base_size = 12)

Interpretasi output: Forest plot IRR menampilkan estimasi titik dan interval kepercayaan 95% untuk setiap koefisien. Titik di kanan garis putus-putus (IRR > 1) berarti prediktor meningkatkan rata-rata hompop. Jika interval kepercayaan tidak melewati garis 1, koefisien signifikan secara statistik.

4.9 Prediksi Rate Poisson

grid_pois <- expand.grid(
  log_inc = seq(min(data_bersih$log_inc),
                max(data_bersih$log_inc),
                length.out = 120),
  race    = factor("Black", levels = levels(data_bersih$race))
)

pred_link <- predict(model_poisson, newdata = grid_pois,
                     type = "link", se.fit = TRUE)

grid_pois_plot <- grid_pois |>
  mutate(
    fit_link  = pred_link$fit,
    se_link   = pred_link$se.fit,
    rate      = exp(fit_link),
    rate_low  = exp(fit_link - 1.96 * se_link),
    rate_high = exp(fit_link + 1.96 * se_link)
  )

ggplot(grid_pois_plot, aes(x = log_inc, y = rate)) +
  geom_ribbon(aes(ymin = rate_low, ymax = rate_high),
              fill = "#2a9d8f", alpha = 0.18) +
  geom_line(color = "#2a9d8f", linewidth = 1.35) +
  labs(title    = "Prediksi Rata-Rata Anggota Rumah Tangga",
       subtitle = "Ras = Black; pita = IK 95%.",
       x        = "Log pendapatan (ln USD)",
       y        = "Prediksi jumlah anggota RT") +
  theme_minimal(base_size = 12)

Interpretasi output: Grafik menampilkan rata-rata prediksi hompop sebagai fungsi log-pendapatan untuk kelompok Black, dengan pita interval kepercayaan 95%. Arah kurva mencerminkan tanda koefisien log_inc.

4.9.1 Prediksi Per Ras

grid_ras <- expand.grid(
  log_inc = seq(min(data_bersih$log_inc),
                max(data_bersih$log_inc),
                length.out = 120),
  race    = factor(c("Black", "White"), levels = levels(data_bersih$race))
)

pred_ras <- predict(model_poisson, newdata = grid_ras,
                    type = "link", se.fit = TRUE)

grid_ras_plot <- grid_ras |>
  mutate(
    rate      = exp(pred_ras$fit),
    rate_low  = exp(pred_ras$fit - 1.96 * pred_ras$se.fit),
    rate_high = exp(pred_ras$fit + 1.96 * pred_ras$se.fit)
  )

ggplot(grid_ras_plot, aes(x = log_inc, y = rate,
                           color = race, fill = race)) +
  geom_ribbon(aes(ymin = rate_low, ymax = rate_high),
              alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.25) +
  scale_color_manual(values = c("Black" = "#e76f51",
                                "White" = "#2a9d8f")) +
  scale_fill_manual(values  = c("Black" = "#e76f51",
                                "White" = "#2a9d8f")) +
  labs(title    = "Prediksi Rata-Rata Anggota RT Menurut Ras",
       subtitle = "Pita = IK 95%.",
       x        = "Log pendapatan (ln USD)",
       y        = "Prediksi jumlah anggota RT",
       color    = "Ras", fill = "Ras") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Interpretasi output: Grafik ini memperlihatkan perbedaan kurva prediksi antara Black dan White. Jarak vertikal antara kedua kurva mencerminkan efek IRR untuk koefisien raceWhite — konstan pada seluruh rentang pendapatan karena model Poisson menggunakan skala log yang additif.

4.10 Pemeriksaan Asumsi Regresi Poisson

Asumsi penting pada regresi Poisson:

1. Respons berupa hitungan nonnegatif — nilai respons harus \(0, 1, 2, \ldots\).
2. Observasi independen — satu observasi tidak bergantung pada observasi lainnya.
3. Mean dimodelkan dengan log link — hubungan prediktor terhadap rata-rata hitungan bersifat linear pada skala log.
4. Equidispersion — secara ideal, \(\text{Var}(Y_i \mid \mathbf{x}_i) = E(Y_i \mid \mathbf{x}_i)\).
5. Exposure ditangani dengan offset bila diperlukan — jika kesempatan terjadinya kejadian berbeda antar-observasi, gunakan offset.

4.10.1 Pemeriksaan Overdispersion

Salah satu ukuran sederhana untuk memeriksa indikasi overdispersion adalah: \[\hat{\phi} = \frac{\sum r_{P,i}^2}{df},\] dengan \(r_{P,i}\) adalah residual Pearson dan \(df\) adalah derajat bebas residual. Nilai \(\hat{\phi}\) sekitar 1 mengindikasikan equidispersion.

phi_hat <- sum(residuals(model_poisson, type = "pearson")^2) /
           df.residual(model_poisson)

disp_tbl <- data.frame(
  `Dispersion Pearson` = round(phi_hat, 3),
  `Interpretasi`       = dplyr::case_when(
    phi_hat < 1.5 ~ "Tidak ada indikasi overdispersion berat",
    phi_hat < 2.5 ~ "Ada indikasi overdispersion sedang",
    TRUE          ~ "Ada indikasi overdispersion kuat — pertimbangkan quasi-Poisson atau NB"
  ),
  check.names = FALSE
)

knitr::kable(disp_tbl,
             caption = "Indikasi overdispersion pada model Poisson")

Indikasi overdispersion pada model Poisson

Dispersion Pearson	Interpretasi
0.838	Tidak ada indikasi overdispersion berat

Interpretasi output: Nilai dispersion Pearson yang mendekati 1 mengindikasikan equidispersion dan model Poisson sudah sesuai. Jika nilai ini jauh lebih besar dari 1, standard error yang dilaporkan akan terlalu kecil, sehingga perlu dipertimbangkan alternatif seperti quasi-Poisson untuk memperbaiki standard error, negative binomial regression ketika varians jauh lebih besar dari mean, atau zero-inflated Poisson ketika jumlah nol jauh lebih banyak dari yang diprediksi.

4.10.2 Plot Residual

data_resid <- data.frame(
  fitted  = fitted(model_poisson),
  pearson = residuals(model_poisson, type = "pearson")
)

ggplot(data_resid, aes(x = fitted, y = pearson)) +
  geom_hline(yintercept = 0, linetype = "dashed",
             color = "#6c757d", linewidth = 0.8) +
  geom_point(alpha = 0.35, size = 1.4, color = "#2a9d8f") +
  geom_smooth(method = "loess", se = FALSE,
              color = "#e76f51", linewidth = 1) +
  labs(title    = "Residual Pearson vs Nilai Fitted",
       subtitle = "Pola acak di sekitar nol mengindikasikan kecocokan model yang baik.",
       x = "Nilai fitted (rata-rata prediksi)",
       y = "Residual Pearson") +
  theme_minimal(base_size = 12)

Interpretasi output: Jika titik-titik tersebar secara acak di sekitar garis nol tanpa pola sistematis, model Poisson sudah cukup sesuai. Pola berbentuk kipas (varians meningkat seiring fitted) mengindikasikan overdispersion.

4.11 Goodness-of-Fit

dev_null <- model_poisson$null.deviance
dev_res  <- model_poisson$deviance
df_null  <- model_poisson$df.null
df_res   <- model_poisson$df.residual
g2       <- dev_null - dev_res
df_g2    <- df_null  - df_res
p_g2     <- 1 - pchisq(g2, df_g2)

gof_tbl <- data.frame(
  Metrik = c("Null deviance", "Residual deviance",
             "G² (penurunan deviance)", "df G²", "p-value G²",
             "AIC", "Pseudo-R² (McFadden)"),
  Nilai  = c(
    round(dev_null, 2),
    round(dev_res,  2),
    round(g2,       2),
    df_g2,
    signif(p_g2, 3),
    round(AIC(model_poisson), 2),
    round(1 - dev_res / dev_null, 4)
  )
)

knitr::kable(gof_tbl, caption = "Ukuran kecocokan model Poisson")

Ukuran kecocokan model Poisson

Metrik	Nilai
Null deviance	9.66860e+02
Residual deviance	9.33720e+02
G² (penurunan deviance)	3.31400e+01
df G²	2.00000e+00
p-value G²	1.00000e-07
AIC	4.26944e+03
Pseudo-R² (McFadden)	3.43000e-02

Interpretasi output: Nilai \(G^2\) menguji apakah model dengan prediktor secara signifikan lebih baik daripada model hanya dengan intercept. p-value yang kecil (< 0,05) berarti prediktor secara bersama-sama berkontribusi signifikan. Pseudo-\(R^2\) McFadden memberikan gambaran proporsi deviance yang berhasil dijelaskan model; nilai 0,2–0,4 sudah dianggap baik untuk model regresi hitung.

4.12 Kesimpulan

Analisis regresi Poisson terhadap data GSS 2024 menghasilkan beberapa temuan utama:

Pertama, formulasi model. Model Poisson membangun satu persamaan log-linear yang menjelaskan rata-rata jumlah anggota rumah tangga (hompop) sebagai fungsi dari log-pendapatan dan ras. Koefisien diinterpretasikan sebagai incidence rate ratio (IRR) setelah dieksponensialkan.

Kedua, peran pendapatan. Log-pendapatan (log_inc) merupakan prediktor numerik utama. Jika koefisiennya negatif dan signifikan, maka rumah tangga dengan pendapatan lebih tinggi cenderung memiliki lebih sedikit anggota — kemungkinan mencerminkan perbedaan gaya hidup atau biaya hidup.

Ketiga, peran ras. Koefisien raceWhite menangkap perbedaan rata-rata ukuran rumah tangga antara kelompok White dan Black setelah mengendalikan pendapatan. Perbedaan ini dapat mencerminkan faktor budaya, historis, maupun struktural yang memengaruhi komposisi rumah tangga.

Keempat, pemeriksaan asumsi. Nilai dispersion Pearson perlu dicermati untuk memastikan asumsi equidispersion terpenuhi. Jika overdispersion terindikasi, quasi-Poisson atau negative binomial regression menjadi alternatif yang lebih robust.

Kelima, keterbatasan. Model ini memberikan gambaran awal yang informatif, namun tidak boleh diinterpretasikan sebagai hubungan kausal. Faktor-faktor lain seperti wilayah geografis, status pernikahan, jumlah kamar rumah, atau tingkat pendidikan dapat memengaruhi ukuran rumah tangga dan belum dimasukkan dalam model ini.