Pengaplikasian Regresi Logistik Multinomial, Regresi Logistik Ordinal, dan Regresi Poisson Menggunakan R

library(dplyr)
library(ggplot2)
library(MASS)
library(nnet)
library(knitr)
library(tidyr)
library(tibble)
library(foreign)
library(scales)

knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.align = "center",
  fig.width = 8.5,
  fig.height = 5,
  dpi = 120,
  out.width = "92%",
  comment = "##"
)

theme_brown <- function(base_size = 12) {
  theme_minimal(base_size = base_size) +
    theme(
      plot.title = element_text(face = "bold", color = "#1f3a4d", size = base_size + 3),
      plot.subtitle = element_text(color = "#65727f", size = base_size),
      axis.title = element_text(color = "#243142", face = "bold"),
      axis.text = element_text(color = "#3f4a55"),
      panel.grid.major = element_line(color = "#e7edf1", linewidth = 0.35),
      panel.grid.minor = element_blank(),
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "white", color = NA),
      plot.margin = margin(14, 16, 12, 16)
    )
}

Pendahuluan

Pengaplikasian Regresi Poisson, Regresi Logistik Ordinal, dan Regresi Logistik Multinomial ini menggunakan software R dan analisis dilakukan menggunakan data provinsi Indonesia dan dataset pendidikan hsbdemo yang juga merupakan data asli dari UCLA Statistical Consulting.

Setiap pembahasan meliputi:

1. Pendahuluan kasus
2. Deskripsi data
3. Eksplorasi data
4. Estimasi model
5. Odds ratio atau rate ratio
6. Prediksi probabilitas atau nilai prediksi
7. Interpretasi hasil
8. Kesimpulan

Sumber Data

Pada ebook ini digunakan tiga studi kasus untuk mengilustrasikan penerapan Regresi Poisson, Regresi Logistik Ordinal, dan Regresi Logistik Multinomial menggunakan perangkat lunak R. Seluruh data yang digunakan berasal dari sumber yang terbuka dan dapat diakses secara publik.

1. Regresi Poisson

Studi kasus Regresi Poisson menggunakan data tingkat provinsi di Indonesia tahun 2024 yang diperoleh dari Badan Pusat Statistik (BPS). Variabel yang digunakan terdiri atas:

• Y (Jumlah Tindak Pidana Menurut Kepolisian Daerah, 2024) https://www.bps.go.id/id/statistics-table/2/MTAxIzI=/jumlah-tindak-pidana-menurut-kepolisian-daerah.html

• X1 (Kepadatan Penduduk, 2024) https://www.bps.go.id/id/statistics-table/3/V1ZSbFRUY3lTbFpEYTNsVWNGcDZjek53YkhsNFFUMDkjMyMwMDAw/jumlah-penduduk--laju-pertumbuhan-penduduk--distribusi-persentase-penduduk--kepadatan-penduduk--rasio-jenis-kelamin-penduduk-menurut-provinsi.html?year=2023

• X2 (Persentase Penduduk Miskin, 2024) https://www.bps.go.id/id/statistics-table/3/UkVkWGJVZFNWakl6VWxKVFQwWjVWeTlSZDNabVFUMDkjMw==/jumlah-dan-persentase-penduduk-miskin-menurut-provinsi--2024.html?year=2024

• X3 (Tingkat Pengangguran Terbuka, 2024) https://www.bps.go.id/id/statistics-table/2/NTQzIzI=/tingkat-pengangguran-terbuka-menurut-provinsi.html

2. Regresi Logistik Ordinal

Studi kasus Regresi Logistik Ordinal menggunakan data tingkat provinsi di Indonesia tahun 2021 yang diperoleh dari Badan Pusat Statistik (BPS). Variabel respon berupa Indeks Kebahagiaan yang kemudian dikategorikan menjadi tiga tingkat, yaitu:

Skor Indeks Kebahagiaan	Kategori
< 70	Rendah
70–75	Sedang
> 75	Tinggi

Variabel yang digunakan terdiri atas:

• Y (Indeks Kebahagiaan Menurut Provinsi, 2021) https://www.bps.go.id/id/statistics-table/2/NjAxIzI=/indeks-kebahagiaan-menurut-provinsi.html

• X1 (Pengeluaran per Kapita Disesuaikan, 2021) https://www.bps.go.id/id/statistics-table/2/NDE2IzI=/-metode-baru-pengeluaran-per-kapita-disesuaikan.html

• X2 (Indeks Pembangunan Manusia (IPM), 2021) https://www.bps.go.id/id/statistics-table/2/NDk0IzI=/-metode-baru-indeks-pembangunan-manusia-menurut-provinsi.html

• X3 (Persentase Penduduk Miskin, 2021) https://www.bps.go.id/id/statistics-table/3/UkVkWGJVZFNWakl6VWxKVFQwWjVWeTlSZDNabVFUMDkjMw==/jumlah-dan-persentase-penduduk-miskin-menurut-provinsi--2024.html?year=2024

3. Regresi Logistik Multinomial

Studi kasus Regresi Logistik Multinomial menggunakan dataset hsbdemo (High School and Beyond) yang tersedia dalam berbagai paket dan sumber pembelajaran statistika.

Sumber data: UCLA Statistical Consulting – hsbdemo Dataset

Variabel yang digunakan terdiri atas:

• Y (Program Pendidikan)
  • General
  • Academic
  • Vocational
• X1 (Nilai Matematika)
• X2 (Nilai Membaca)
• X3 (Nilai Menulis)
• X4 (Status Sosial Ekonomi)

Dataset ini sering digunakan sebagai contoh pembelajaran regresi logistik multinomial karena memiliki variabel respon kategorik dengan lebih dari dua kategori yang tidak memiliki urutan tertentu.

Cek Struktur Data Asli

Data provinsi dibaca dari file CSV asli. File ini memakai separator titik koma dan koma sebagai tanda desimal, sehingga fungsi yang sesuai adalah read.csv2().

file_data <- "D:/Semester 4/ADK/Data ADK 1_Regresi poisson dan regresi ordinal.csv"

stopifnot(file.exists(file_data))

data <- read.csv2(
  file = file_data,
  header = TRUE,
  check.names = FALSE,
  stringsAsFactors = FALSE
)

names(data)[1] <- "Provinsi"

dim(data)

## [1] 35  9

names(data)

## [1] "Provinsi"                                                    
## [2] "Jumlah Tindak Pidana Menurut Kepolisian Daerah"              
## [3] "Kepadatan Penduduk per km persegi (Km2)"                     
## [4] "Persentase Penduduk Miskin (2024)"                           
## [5] "Tingkat Pengangguran Terbuka"                                
## [6] "Indeks Kebahagiaan Menurut Provinsi (2021)"                  
## [7] "Pengeluaran per Kapita Disesuaikan (Ribu Rupiah/Orang/Tahun)"
## [8] "Persentase Penduduk Miskin (2021)"                           
## [9] "Indeks Pembangunan Manusia (IPM) menurut Provinsi"

str(data)

## 'data.frame':    35 obs. of  9 variables:
##  $ Provinsi                                                    : chr  "ACEH" "SUMATERA UTARA" "SUMATERA BARAT" "RIAU" ...
##  $ Jumlah Tindak Pidana Menurut Kepolisian Daerah              : int  11934 60724 13099 17799 7592 25154 5323 16249 2645 5294 ...
##  $ Kepadatan Penduduk per km persegi (Km2)                     : int  98 215 139 75 76 102 105 281 92 264 ...
##  $ Persentase Penduduk Miskin (2024)                           : num  12.64 7.19 5.42 6.36 7.26 ...
##  $ Tingkat Pengangguran Terbuka                                : num  5.75 5.6 5.75 3.7 4.48 3.86 3.11 4.19 4.63 6.39 ...
##  $ Indeks Kebahagiaan Menurut Provinsi (2021)                  : chr  "Sedang" "Sedang" "Sedang" "Sedang" ...
##  $ Pengeluaran per Kapita Disesuaikan (Ribu Rupiah/Orang/Tahun): int  9572 10499 10790 10736 10588 10662 10487 10038 12819 14122 ...
##  $ Persentase Penduduk Miskin (2021)                           : num  15.33 9.01 6.63 7.12 8.09 ...
##  $ Indeks Pembangunan Manusia (IPM) menurut Provinsi           : num  72.2 72 72.7 72.9 71.6 ...

missing_value <- colSums(is.na(data))
missing_value

##                                                     Provinsi 
##                                                            0 
##               Jumlah Tindak Pidana Menurut Kepolisian Daerah 
##                                                            0 
##                      Kepadatan Penduduk per km persegi (Km2) 
##                                                            0 
##                            Persentase Penduduk Miskin (2024) 
##                                                            0 
##                                 Tingkat Pengangguran Terbuka 
##                                                            0 
##                   Indeks Kebahagiaan Menurut Provinsi (2021) 
##                                                            0 
## Pengeluaran per Kapita Disesuaikan (Ribu Rupiah/Orang/Tahun) 
##                                                            0 
##                            Persentase Penduduk Miskin (2021) 
##                                                            0 
##            Indeks Pembangunan Manusia (IPM) menurut Provinsi 
##                                                            0

data %>%
  dplyr::count(`Indeks Kebahagiaan Menurut Provinsi (2021)`) %>%
  kable(caption = "Distribusi kategori indeks kebahagiaan")

Distribusi kategori indeks kebahagiaan

Indeks Kebahagiaan Menurut Provinsi (2021)	n
Rendah	4
Sedang	27
Tinggi	4

Model 1

Regresi Poisson

BAB 1 REGRESI POISSON

1.1 Pendahuluan Kasus

Kasus pada penelitian ini membahas faktor-faktor yang mempengaruhi jumlah tindak pidana menurut kepolisian daerah di Indonesia tahun 2024. Variabel respon yang digunakan adalah jumlah tindak pidana, sedangkan variabel prediktor meliputi kepadatan penduduk, persentase penduduk miskin, dan tingkat pengangguran terbuka.

Karena variabel respon berbentuk count data atau data cacahan, maka metode yang digunakan adalah Regresi Poisson.

1.2 Deskripsi Data

poisson_data <- data %>%
  dplyr::select(
    Provinsi,
    Tindak_Pidana = `Jumlah Tindak Pidana Menurut Kepolisian Daerah`,
    Kepadatan_Penduduk = `Kepadatan Penduduk per km persegi (Km2)`,
    Persentase_Miskin = `Persentase Penduduk Miskin (2024)`,
    TPT = `Tingkat Pengangguran Terbuka`
  )

head(poisson_data)

##           Provinsi Tindak_Pidana Kepadatan_Penduduk Persentase_Miskin  TPT
## 1             ACEH         11934                 98             12.64 5.75
## 2   SUMATERA UTARA         60724                215              7.19 5.60
## 3   SUMATERA BARAT         13099                139              5.42 5.75
## 4             RIAU         17799                 75              6.36 3.70
## 5            JAMBI          7592                 76              7.26 4.48
## 6 SUMATERA SELATAN         25154                102             10.51 3.86

str(poisson_data)

## 'data.frame':    35 obs. of  5 variables:
##  $ Provinsi          : chr  "ACEH" "SUMATERA UTARA" "SUMATERA BARAT" "RIAU" ...
##  $ Tindak_Pidana     : int  11934 60724 13099 17799 7592 25154 5323 16249 2645 5294 ...
##  $ Kepadatan_Penduduk: int  98 215 139 75 76 102 105 281 92 264 ...
##  $ Persentase_Miskin : num  12.64 7.19 5.42 6.36 7.26 ...
##  $ TPT               : num  5.75 5.6 5.75 3.7 4.48 3.86 3.11 4.19 4.63 6.39 ...

summary(poisson_data)

##    Provinsi         Tindak_Pidana    Kepadatan_Penduduk Persentase_Miskin
##  Length:35          Min.   :  1593   Min.   :   10.0    Min.   : 3.800   
##  Class :character   1st Qu.:  5631   1st Qu.:   63.0    1st Qu.: 5.605   
##  Mode  :character   Median :  9623   Median :  105.0    Median : 7.770   
##                     Mean   : 32114   Mean   :  738.5    Mean   : 9.151   
##                     3rd Qu.: 17024   3rd Qu.:  272.5    3rd Qu.:10.875   
##                     Max.   :561993   Max.   :16165.0    Max.   :21.090   
##       TPT       
##  Min.   :1.790  
##  1st Qu.:3.590  
##  Median :4.190  
##  Mean   :4.478  
##  3rd Qu.:5.675  
##  Max.   :6.750

poisson_data %>%
  dplyr::summarise(
    n = dplyr::n(),
    rata_rata_tindak_pidana = mean(Tindak_Pidana),
    median_tindak_pidana = median(Tindak_Pidana),
    minimum = min(Tindak_Pidana),
    maksimum = max(Tindak_Pidana)
  ) %>%
  kable(digits = 3, caption = "Ringkasan jumlah tindak pidana")

Ringkasan jumlah tindak pidana

n	rata_rata_tindak_pidana	median_tindak_pidana	minimum	maksimum
35	32113.89	9623	1593	561993

1.3 Eksplorasi Data

ggplot(
  poisson_data,
  aes(x = Kepadatan_Penduduk, y = Tindak_Pidana)
) +
  geom_point(color = "#2f7f73", size = 2.6, alpha = 0.9) +
  geom_smooth(method = "lm", se = FALSE, color = "#d18b2f") +
  labs(
    title = "Jumlah Tindak Pidana dan Kepadatan Penduduk",
    subtitle = "Setiap titik merepresentasikan satu provinsi.",
    x = "Kepadatan penduduk per km persegi",
    y = "Jumlah tindak pidana"
  ) +
  theme_brown()

cor(poisson_data[, 2:5]) %>%
  round(3) %>%
  kable(caption = "Matriks korelasi variabel model Poisson")

Matriks korelasi variabel model Poisson

	Tindak_Pidana	Kepadatan_Penduduk	Persentase_Miskin	TPT
Tindak_Pidana	1.000	0.084	-0.053	0.113
Kepadatan_Penduduk	0.084	1.000	-0.222	0.256
Persentase_Miskin	-0.053	-0.222	1.000	-0.163
TPT	0.113	0.256	-0.163	1.000

1.4 Estimasi Model

model_poisson <- glm(
  Tindak_Pidana ~ Kepadatan_Penduduk + Persentase_Miskin + TPT,
  family = poisson(link = "log"),
  data = poisson_data
)

summary(model_poisson)

## 
## Call:
## glm(formula = Tindak_Pidana ~ Kepadatan_Penduduk + Persentase_Miskin + 
##     TPT, family = poisson(link = "log"), data = poisson_data)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         9.485e+00  4.571e-03 2075.26   <2e-16 ***
## Kepadatan_Penduduk  2.813e-05  2.572e-07  109.37   <2e-16 ***
## Persentase_Miskin  -1.959e-02  2.383e-04  -82.21   <2e-16 ***
## TPT                 2.210e-01  7.808e-04  283.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3006885  on 34  degrees of freedom
## Residual deviance: 2864818  on 31  degrees of freedom
## AIC: 2865218
## 
## Number of Fisher Scoring iterations: 7

pois_coef <- as.data.frame(coef(summary(model_poisson))) %>%
  tibble::rownames_to_column("parameter") %>%
  dplyr::rename(
    estimate = Estimate,
    std_error = `Std. Error`,
    z_value = `z value`,
    p_value = `Pr(>|z|)`
  ) %>%
  dplyr::mutate(
    IRR = exp(estimate),
    CI_low = exp(estimate - 1.96 * std_error),
    CI_high = exp(estimate + 1.96 * std_error)
  )

pois_coef %>%
  dplyr::mutate(
    dplyr::across(
      c(estimate, std_error, z_value, p_value, IRR, CI_low, CI_high),
      ~ round(.x, 4)
    )
  ) %>%
  kable(
    caption = "Ringkasan hasil regresi Poisson",
    col.names = c(
      "Parameter", "Estimate", "SE", "z-value", "p-value",
      "IRR", "CI 2.5%", "CI 97.5%"
    )
  )

Ringkasan hasil regresi Poisson

Parameter	Estimate	SE	z-value	p-value	IRR	CI 2.5%	CI 97.5%
(Intercept)	9.4853	0.0046	2075.2650	0	13164.5342	13047.1271	13282.9978
Kepadatan_Penduduk	0.0000	0.0000	109.3689	0	1.0000	1.0000	1.0000
Persentase_Miskin	-0.0196	0.0002	-82.2081	0	0.9806	0.9801	0.9811
TPT	0.2210	0.0008	283.1217	0	1.2474	1.2455	1.2493

1.5 Rate Ratio

exp(coef(model_poisson)) %>%
  round(4)

##        (Intercept) Kepadatan_Penduduk  Persentase_Miskin                TPT 
##         13164.5342             1.0000             0.9806             1.2474

1.6 Prediksi Nilai Respons

poisson_data <- poisson_data %>%
  dplyr::mutate(
    Prediksi = predict(model_poisson, type = "response"),
    Residual = Tindak_Pidana - Prediksi
  )

head(poisson_data)

##           Provinsi Tindak_Pidana Kepadatan_Penduduk Persentase_Miskin  TPT
## 1             ACEH         11934                 98             12.64 5.75
## 2   SUMATERA UTARA         60724                215              7.19 5.60
## 3   SUMATERA BARAT         13099                139              5.42 5.75
## 4             RIAU         17799                 75              6.36 3.70
## 5            JAMBI          7592                 76              7.26 4.48
## 6 SUMATERA SELATAN         25154                102             10.51 3.86
##   Prediksi     Residual
## 1 36734.39 -24800.38611
## 2 39670.48  21053.51501
## 3 42364.02 -29265.02392
## 4 26388.64  -8589.64080
## 5 30807.26 -23215.25946
## 6 25223.26    -69.26116

ggplot(poisson_data, aes(x = Tindak_Pidana, y = Prediksi)) +
  geom_point(color = "#5568b8", size = 2.6, alpha = 0.9) +
  geom_abline(slope = 1, intercept = 0, color = "#b84a3a", linewidth = 0.8) +
  labs(
    title = "Aktual dan Prediksi Model Poisson",
    subtitle = "Garis merah menunjukkan prediksi yang sama dengan nilai aktual.",
    x = "Jumlah tindak pidana aktual",
    y = "Jumlah tindak pidana prediksi"
  ) +
  theme_brown()

1.7 Interpretasi Hasil

Berdasarkan model Poisson, koefisien dapat dibaca melalui nilai IRR = exp(beta). Nilai IRR di atas 1 menunjukkan peningkatan rata-rata jumlah tindak pidana, sedangkan nilai IRR di bawah 1 menunjukkan penurunan rata-rata jumlah tindak pidana, dengan prediktor lain konstan.

Pada data ini, kepadatan penduduk dan tingkat pengangguran terbuka memiliki arah hubungan positif terhadap rata-rata jumlah tindak pidana. Persentase penduduk miskin memiliki arah hubungan negatif terhadap rata-rata jumlah tindak pidana.

1.8 Kesimpulan

Regresi Poisson mampu digunakan untuk menganalisis faktor-faktor yang mempengaruhi jumlah tindak pidana di Indonesia. Model ini sesuai karena variabel respon berbentuk hitungan.

Model 2

Regresi Logistik Ordinal

BAB 2 REGRESI LOGISTIK ORDINAL

2.1 Pendahuluan Kasus

Kasus ini membahas faktor-faktor yang mempengaruhi indeks kebahagiaan menurut provinsi di Indonesia tahun 2021. Variabel respon berbentuk kategori ordinal yaitu rendah, sedang, dan tinggi.

Karena variabel respon berbentuk ordinal, maka digunakan Regresi Logistik Ordinal.

2.2 Deskripsi Data

ordinal_data <- data %>%
  dplyr::select(
    Provinsi,
    Kebahagiaan = `Indeks Kebahagiaan Menurut Provinsi (2021)`,
    Pengeluaran = `Pengeluaran per Kapita Disesuaikan (Ribu Rupiah/Orang/Tahun)`,
    Miskin = `Persentase Penduduk Miskin (2021)`,
    IPM = `Indeks Pembangunan Manusia (IPM) menurut Provinsi`
  ) %>%
  dplyr::mutate(
    Kebahagiaan = factor(
      Kebahagiaan,
      levels = c("Rendah", "Sedang", "Tinggi"),
      ordered = TRUE
    ),
    Pengeluaran_z = as.numeric(scale(Pengeluaran)),
    Miskin_z = as.numeric(scale(Miskin)),
    IPM_z = as.numeric(scale(IPM))
  )

head(ordinal_data)

##           Provinsi Kebahagiaan Pengeluaran Miskin   IPM Pengeluaran_z
## 1             ACEH      Sedang        9572  15.33 72.18  -0.549371808
## 2   SUMATERA UTARA      Sedang       10499   9.01 72.00  -0.124557272
## 3   SUMATERA BARAT      Sedang       10790   6.63 72.65   0.008798748
## 4             RIAU      Sedang       10736   7.12 72.94  -0.015947730
## 5            JAMBI      Tinggi       10588   8.09 71.63  -0.083771410
## 6 SUMATERA SELATAN      Sedang       10662  12.84 70.24  -0.049859570
##     Miskin_z       IPM_z
## 1  0.8613277  0.20407745
## 2 -0.3258222  0.15771294
## 3 -0.7728818  0.32514034
## 4 -0.6808401  0.39983873
## 5 -0.4986352  0.06240811
## 6  0.3936057 -0.29562896

summary(ordinal_data)

##    Provinsi         Kebahagiaan  Pengeluaran        Miskin      
##  Length:35          Rendah: 4   Min.   : 6955   Min.   : 4.530  
##  Class :character   Sedang:27   1st Qu.: 9380   1st Qu.: 6.775  
##  Mode  :character   Tinggi: 4   Median :10662   Median : 9.010  
##                                 Mean   :10771   Mean   :10.745  
##                                 3rd Qu.:11446   3rd Qu.:12.920  
##                                 Max.   :18520   Max.   :26.860  
##       IPM        Pengeluaran_z         Miskin_z           IPM_z         
##  Min.   :60.62   Min.   :-1.74866   Min.   :-1.1673   Min.   :-2.77355  
##  1st Qu.:69.75   1st Qu.:-0.63759   1st Qu.:-0.7456   1st Qu.:-0.42184  
##  Median :71.66   Median :-0.04986   Median :-0.3258   Median : 0.07014  
##  Mean   :71.39   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.:72.55   3rd Qu.: 0.30919   3rd Qu.: 0.4086   3rd Qu.: 0.29938  
##  Max.   :81.11   Max.   : 3.55121   Max.   : 3.0271   Max.   : 2.50427

2.3 Eksplorasi Data

ordinal_data %>%
  dplyr::count(Kebahagiaan) %>%
  dplyr::mutate(proporsi = n / sum(n)) %>%
  kable(digits = 3, caption = "Distribusi kategori kebahagiaan")

Distribusi kategori kebahagiaan

Kebahagiaan	n	proporsi
Rendah	4	0.114
Sedang	27	0.771
Tinggi	4	0.114

ggplot(ordinal_data, aes(x = Kebahagiaan, y = Pengeluaran, fill = Kebahagiaan)) +
  geom_boxplot(width = 0.65, alpha = 0.88) +
  scale_fill_manual(values = c("#b84a3a", "#d18b2f", "#2f7f73")) +
  labs(
    title = "Pengeluaran per Kapita Menurut Kategori Kebahagiaan",
    subtitle = "Respons ordinal memiliki urutan Rendah, Sedang, dan Tinggi.",
    x = "Kategori kebahagiaan",
    y = "Pengeluaran per kapita disesuaikan"
  ) +
  theme_brown() +
  theme(legend.position = "none")

2.4 Estimasi Model

Prediktor numerik pada model ordinal distandarkan agar proses estimasi lebih stabil. Pendekatan analisis tetap sama, yaitu model proportional odds dengan fungsi MASS::polr().

model_ordinal <- MASS::polr(
  Kebahagiaan ~ Pengeluaran_z + Miskin_z + IPM_z,
  data = ordinal_data,
  Hess = TRUE
)

summary(model_ordinal)

## Call:
## MASS::polr(formula = Kebahagiaan ~ Pengeluaran_z + Miskin_z + 
##     IPM_z, data = ordinal_data, Hess = TRUE)
## 
## Coefficients:
##                Value Std. Error t value
## Pengeluaran_z -2.566     1.0429  -2.461
## Miskin_z      -1.089     0.6284  -1.733
## IPM_z          1.899     1.0779   1.762
## 
## Intercepts:
##               Value   Std. Error t value
## Rendah|Sedang -2.7369  0.7251    -3.7746
## Sedang|Tinggi  2.6697  0.6854     3.8951
## 
## Residual Deviance: 38.44048 
## AIC: 48.44048

ordinal_coef <- as.data.frame(coef(summary(model_ordinal))) %>%
  tibble::rownames_to_column("parameter") %>%
  dplyr::rename(
    estimate = Value,
    std_error = `Std. Error`,
    t_value = `t value`
  ) %>%
  dplyr::mutate(
    p_value = 2 * (1 - pnorm(abs(t_value))),
    OR = exp(estimate),
    CI_low = exp(estimate - 1.96 * std_error),
    CI_high = exp(estimate + 1.96 * std_error)
  )

ordinal_coef %>%
  dplyr::mutate(
    dplyr::across(
      c(estimate, std_error, t_value, p_value, OR, CI_low, CI_high),
      ~ round(.x, 4)
    )
  ) %>%
  kable(
    caption = "Ringkasan hasil regresi logistik ordinal",
    col.names = c(
      "Parameter", "Estimate", "SE", "t-value", "p-value",
      "OR", "CI 2.5%", "CI 97.5%"
    )
  )

Ringkasan hasil regresi logistik ordinal

Parameter	Estimate	SE	t-value	p-value	OR	CI 2.5%	CI 97.5%
Pengeluaran_z	-2.5664	1.0429	-2.4608	0.0139	0.0768	0.0099	0.5931
Miskin_z	-1.0894	0.6284	-1.7334	0.0830	0.3364	0.0982	1.1530
IPM_z	1.8988	1.0779	1.7616	0.0781	6.6776	0.8075	55.2212
Rendah|Sedang	-2.7369	0.7251	-3.7746	0.0002	0.0648	0.0156	0.2683
Sedang|Tinggi	2.6697	0.6854	3.8951	0.0001	14.4361	3.7672	55.3203

2.5 Odds Ratio

exp(coef(model_ordinal)) %>%
  round(4)

## Pengeluaran_z      Miskin_z         IPM_z 
##        0.0768        0.3364        6.6776

2.6 Prediksi Probabilitas

prediksi_ordinal <- predict(model_ordinal, type = "probs")

head(prediksi_ordinal)

##       Rendah    Sedang     Tinggi
## 1 0.02669959 0.8327389 0.14056149
## 2 0.02386935 0.8211000 0.15503066
## 3 0.01516224 0.7591827 0.22565507
## 4 0.01367057 0.7417867 0.24454271
## 5 0.02624655 0.8310549 0.14269858
## 6 0.13298738 0.8385940 0.02841865

ordinal_prediksi <- ordinal_data %>%
  dplyr::select(Provinsi, Kebahagiaan) %>%
  dplyr::bind_cols(as.data.frame(prediksi_ordinal)) %>%
  dplyr::mutate(Kelas_Prediksi = predict(model_ordinal, type = "class"))

head(ordinal_prediksi)

##           Provinsi Kebahagiaan     Rendah    Sedang     Tinggi Kelas_Prediksi
## 1             ACEH      Sedang 0.02669959 0.8327389 0.14056149         Sedang
## 2   SUMATERA UTARA      Sedang 0.02386935 0.8211000 0.15503066         Sedang
## 3   SUMATERA BARAT      Sedang 0.01516224 0.7591827 0.22565507         Sedang
## 4             RIAU      Sedang 0.01367057 0.7417867 0.24454271         Sedang
## 5            JAMBI      Tinggi 0.02624655 0.8310549 0.14269858         Sedang
## 6 SUMATERA SELATAN      Sedang 0.13298738 0.8385940 0.02841865         Sedang

2.7 Interpretasi Hasil

Pada output MASS::polr(), tanda koefisien mengikuti konvensi model proportional odds yang digunakan oleh fungsi tersebut. Nilai odds ratio di atas 1 menunjukkan peningkatan odds berada pada kategori kebahagiaan yang lebih tinggi, sedangkan nilai di bawah 1 menunjukkan penurunan odds, untuk setiap kenaikan satu simpangan baku prediktor.

2.8 Kesimpulan

Regresi logistik ordinal dapat digunakan karena variabel respon mempunyai urutan kategori. Pada data ini, model menghasilkan probabilitas prediksi untuk kategori Rendah, Sedang, dan Tinggi pada setiap provinsi.

Model 3

Regresi Logistik Multinomial

BAB 3 REGRESI LOGISTIK MULTINOMIAL

3.1 Pendahuluan Kasus

Kasus ini membahas faktor-faktor yang mempengaruhi pemilihan program pendidikan siswa menggunakan dataset hsbdemo.

Variabel respon terdiri atas tiga kategori:

1. general
2. academic
3. vocation

Karena variabel respon memiliki lebih dari dua kategori nominal, maka digunakan Regresi Logistik Multinomial.

3.2 Deskripsi Data

file_hsbdemo <- "C:/Users/dell/Documents/Codex/2026-06-02/files-mentioned-by-the-user-ebook/outputs/hsbdemo.csv"
url_hsbdemo <- "https://stats.idre.ucla.edu/stat/data/hsbdemo.dta"

if (file.exists(file_hsbdemo)) {
  hsbdemo <- read.csv(file_hsbdemo, stringsAsFactors = TRUE)
} else {
  temp_hsbdemo <- tempfile(fileext = ".dta")
  download.file(url_hsbdemo, temp_hsbdemo, mode = "wb", quiet = TRUE)
  hsbdemo <- foreign::read.dta(temp_hsbdemo)
}

hsbdemo <- hsbdemo %>%
  dplyr::mutate(
    prog = factor(prog, levels = c("general", "academic", "vocation")),
    ses = factor(ses, levels = c("low", "middle", "high"))
  )

multinom_data <- hsbdemo %>%
  dplyr::mutate(
    prog = relevel(prog, ref = "general"),
    ses = relevel(ses, ref = "low")
  )

head(multinom_data)

##    id female    ses schtyp     prog read write math science socst       honors
## 1  45 female    low public vocation   34    35   41      29    26 not enrolled
## 2 108   male middle public  general   34    33   41      36    36 not enrolled
## 3  15   male   high public vocation   39    39   44      26    42 not enrolled
## 4  67   male    low public vocation   37    37   42      33    32 not enrolled
## 5 153   male middle public vocation   39    31   40      39    51 not enrolled
## 6  51 female   high public  general   42    36   42      31    39 not enrolled
##   awards cid
## 1      0   1
## 2      0   1
## 3      0   1
## 4      0   1
## 5      0   1
## 6      0   1

str(multinom_data)

## 'data.frame':    200 obs. of  13 variables:
##  $ id     : int  45 108 15 67 153 51 164 133 2 53 ...
##  $ female : Factor w/ 2 levels "female","male": 1 2 2 2 2 1 2 2 1 2 ...
##  $ ses    : Factor w/ 3 levels "low","middle",..: 1 2 3 1 2 3 2 2 2 2 ...
##  $ schtyp : Factor w/ 2 levels "private","public": 2 2 2 2 2 2 2 2 2 2 ...
##  $ prog   : Factor w/ 3 levels "general","academic",..: 3 1 3 3 3 1 3 3 3 3 ...
##  $ read   : int  34 34 39 37 39 42 31 50 39 34 ...
##  $ write  : int  35 33 39 37 31 36 36 31 41 37 ...
##  $ math   : int  41 41 44 42 40 42 46 40 33 46 ...
##  $ science: int  29 36 26 33 39 31 39 34 42 39 ...
##  $ socst  : int  26 36 42 32 51 39 46 31 41 31 ...
##  $ honors : Factor w/ 2 levels "enrolled","not enrolled": 2 2 2 2 2 2 2 2 2 2 ...
##  $ awards : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ cid    : int  1 1 1 1 1 1 1 1 1 1 ...

summary(multinom_data)

##        id            female        ses         schtyp          prog    
##  Min.   :  1.00   female:109   low   :47   private: 32   general : 45  
##  1st Qu.: 50.75   male  : 91   middle:95   public :168   academic:105  
##  Median :100.50                high  :58                 vocation: 50  
##  Mean   :100.50                                                        
##  3rd Qu.:150.25                                                        
##  Max.   :200.00                                                        
##       read           write            math          science     
##  Min.   :28.00   Min.   :31.00   Min.   :33.00   Min.   :26.00  
##  1st Qu.:44.00   1st Qu.:45.75   1st Qu.:45.00   1st Qu.:44.00  
##  Median :50.00   Median :54.00   Median :52.00   Median :53.00  
##  Mean   :52.23   Mean   :52.77   Mean   :52.65   Mean   :51.85  
##  3rd Qu.:60.00   3rd Qu.:60.00   3rd Qu.:59.00   3rd Qu.:58.00  
##  Max.   :76.00   Max.   :67.00   Max.   :75.00   Max.   :74.00  
##      socst                honors        awards          cid       
##  Min.   :26.00   enrolled    : 53   Min.   :0.00   Min.   : 1.00  
##  1st Qu.:46.00   not enrolled:147   1st Qu.:0.00   1st Qu.: 5.00  
##  Median :52.00                      Median :1.00   Median :10.50  
##  Mean   :52.41                      Mean   :1.67   Mean   :10.43  
##  3rd Qu.:61.00                      3rd Qu.:2.00   3rd Qu.:15.00  
##  Max.   :71.00                      Max.   :7.00   Max.   :20.00

3.3 Eksplorasi Data

multinom_data %>%
  dplyr::count(prog) %>%
  dplyr::mutate(proporsi = n / sum(n)) %>%
  kable(digits = 3, caption = "Distribusi program pendidikan")

Distribusi program pendidikan

prog	n	proporsi
general	45	0.225
academic	105	0.525
vocation	50	0.250

ggplot(multinom_data, aes(x = prog, y = math, fill = prog)) +
  geom_boxplot(width = 0.65, alpha = 0.88) +
  scale_fill_manual(values = c("#2f7f73", "#5568b8", "#d18b2f")) +
  labs(
    title = "Nilai Matematika Menurut Program Pendidikan",
    subtitle = "Respons multinomial tidak memiliki urutan alami.",
    x = "Program pendidikan",
    y = "Nilai matematika"
  ) +
  theme_brown() +
  theme(legend.position = "none")

3.4 Estimasi Model

model_multinom <- nnet::multinom(
  prog ~ math + read + write + ses,
  data = multinom_data,
  trace = FALSE
)

summary(model_multinom)

## Call:
## nnet::multinom(formula = prog ~ math + read + write + ses, data = multinom_data, 
##     trace = FALSE)
## 
## Coefficients:
##          (Intercept)        math        read       write sesmiddle   seshigh
## academic   -4.590817  0.05982951  0.02274116  0.01207661 0.3245517 0.8622315
## vocation    4.058905 -0.04062592 -0.01856355 -0.03454104 0.9746126 0.3709773
## 
## Std. Errors:
##          (Intercept)       math       read      write sesmiddle   seshigh
## academic    1.390563 0.03007151 0.02647349 0.02669394 0.4615532 0.5370431
## vocation    1.580482 0.03542358 0.03077524 0.02866507 0.5039758 0.6643516
## 
## Residual Deviance: 339.5503 
## AIC: 363.5503

multi_sum <- summary(model_multinom)

coef_multi <- as.data.frame(multi_sum$coefficients)
se_multi <- as.data.frame(multi_sum$standard.errors)

coef_long <- coef_multi %>%
  tibble::rownames_to_column("kategori") %>%
  tidyr::pivot_longer(
    cols = -kategori,
    names_to = "variabel",
    values_to = "estimate"
  )

se_long <- se_multi %>%
  tibble::rownames_to_column("kategori") %>%
  tidyr::pivot_longer(
    cols = -kategori,
    names_to = "variabel",
    values_to = "std_error"
  )

result_multi <- coef_long %>%
  dplyr::left_join(se_long, by = c("kategori", "variabel")) %>%
  dplyr::mutate(
    z_value = estimate / std_error,
    p_value = 2 * (1 - pnorm(abs(z_value))),
    RRR = exp(estimate),
    CI_low = exp(estimate - 1.96 * std_error),
    CI_high = exp(estimate + 1.96 * std_error)
  )

result_multi %>%
  dplyr::mutate(
    dplyr::across(
      c(estimate, std_error, z_value, p_value, RRR, CI_low, CI_high),
      ~ round(.x, 4)
    )
  ) %>%
  kable(
    caption = "Ringkasan hasil regresi logistik multinomial",
    col.names = c(
      "Kategori", "Variabel", "Estimate", "SE", "z-value",
      "p-value", "RRR", "CI 2.5%", "CI 97.5%"
    )
  )

Ringkasan hasil regresi logistik multinomial

Kategori	Variabel	Estimate	SE	z-value	p-value	RRR	CI 2.5%	CI 97.5%
academic	(Intercept)	-4.5908	1.3906	-3.3014	0.0010	0.0101	0.0007	0.1548
academic	math	0.0598	0.0301	1.9896	0.0466	1.0617	1.0009	1.1261
academic	read	0.0227	0.0265	0.8590	0.3903	1.0230	0.9713	1.0775
academic	write	0.0121	0.0267	0.4524	0.6510	1.0121	0.9606	1.0665
academic	sesmiddle	0.3246	0.4616	0.7032	0.4819	1.3834	0.5598	3.4185
academic	seshigh	0.8622	0.5370	1.6055	0.1084	2.3684	0.8267	6.7858
vocation	(Intercept)	4.0589	1.5805	2.5681	0.0102	57.9109	2.6147	1282.6063
vocation	math	-0.0406	0.0354	-1.1469	0.2514	0.9602	0.8958	1.0292
vocation	read	-0.0186	0.0308	-0.6032	0.5464	0.9816	0.9241	1.0426
vocation	write	-0.0345	0.0287	-1.2050	0.2282	0.9660	0.9133	1.0219
vocation	sesmiddle	0.9746	0.5040	1.9338	0.0531	2.6501	0.9869	7.1164
vocation	seshigh	0.3710	0.6644	0.5584	0.5766	1.4492	0.3941	5.3287

3.5 Odds Ratio

exp(coef(model_multinom)) %>%
  round(4)

##          (Intercept)   math   read  write sesmiddle seshigh
## academic      0.0101 1.0617 1.0230 1.0121    1.3834  2.3684
## vocation     57.9109 0.9602 0.9816 0.9660    2.6501  1.4492

3.6 Prediksi Probabilitas

prediksi_multinom <- predict(model_multinom, type = "probs")

head(prediksi_multinom)

##     general   academic  vocation
## 1 0.3196329 0.12461726 0.5557498
## 2 0.1547056 0.08145055 0.7638439
## 3 0.2457138 0.31924970 0.4350365
## 4 0.3415735 0.15506477 0.5033618
## 5 0.1523490 0.08262971 0.7650213
## 6 0.2378192 0.28305785 0.4791230

multinom_prediksi <- multinom_data %>%
  dplyr::select(id, prog, math, read, write, ses) %>%
  dplyr::bind_cols(as.data.frame(prediksi_multinom)) %>%
  dplyr::mutate(Kelas_Prediksi = predict(model_multinom, type = "class"))

head(multinom_prediksi)

##    id     prog math read write    ses   general   academic  vocation
## 1  45 vocation   41   34    35    low 0.3196329 0.12461726 0.5557498
## 2 108  general   41   34    33 middle 0.1547056 0.08145055 0.7638439
## 3  15 vocation   44   39    39   high 0.2457138 0.31924970 0.4350365
## 4  67 vocation   42   37    37    low 0.3415735 0.15506477 0.5033618
## 5 153 vocation   40   39    31 middle 0.1523490 0.08262971 0.7650213
## 6  51  general   42   42    36   high 0.2378192 0.28305785 0.4791230
##   Kelas_Prediksi
## 1       vocation
## 2       vocation
## 3       vocation
## 4       vocation
## 5       vocation
## 6       vocation

3.7 Interpretasi Hasil

Pada regresi logistik multinomial, kategori referensi untuk prog adalah general. Nilai RRR = exp(beta) menunjukkan perubahan relative risk ratio untuk memilih kategori tertentu dibandingkan kategori referensi, dengan prediktor lain konstan.

3.8 Kesimpulan

Regresi logistik multinomial mampu digunakan untuk menganalisis faktor-faktor yang mempengaruhi pemilihan program pendidikan siswa dengan kategori respon lebih dari dua dan tidak berurutan.

Penutup

Secara umum, metode regresi yang digunakan harus disesuaikan dengan bentuk variabel respon.

1. Regresi Poisson digunakan untuk data count.
2. Regresi Logistik Ordinal digunakan untuk data kategori bertingkat.
3. Regresi Logistik Multinomial digunakan untuk data kategori nominal lebih dari dua.

Penggunaan software R mempermudah proses pengolahan data, estimasi model, interpretasi hasil, dan visualisasi data.

Daftar Pustaka

1. Agresti, A. (2018). Statistical Methods for the Social Sciences.
2. Hosmer, D. W. & Lemeshow, S. (2013). Applied Logistic Regression.
3. Badan Pusat Statistik. https://www.bps.go.id
4. UCLA Statistical Consulting. hsbdemo dataset.