Multinomial Logistic Regression pada Data ISPU DKI Jakarta 2021

Generalized Linear Models

Multinomial Logistic Regression

1. Persiapan Data

# Load dataset
df_MAR <- read.csv("C:/Users/Alfathrindra Agastyo/Documents/Semester 4/Analisis Multivariat/Tugas/ISPU DKI JAKARTA 2021.csv")

# Lihat struktur data
View(df_MAR)

Pada tahap ini, dataset yang berisi informasi tentang indeks standar pencemar udara (ISPU) di DKI Jakarta tahun 2021 dibaca ke dalam R dan ditampilkan strukturnya untuk memahami tipe data dari masing-masing kolom.

2. Transformasi Variabel Kategori

library(tidyverse)

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

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

# Transformasi variabel kategori menjadi faktor dengan label yang lebih deskriptif
df_MAR <- df_MAR %>%
  mutate(
    categori = factor(categori, 
                      labels = c("SEDANG", "TIDAK SEHAT", "BAIK"))
  )

# Periksa struktur data setelah transformasi
str(df_MAR)

## 'data.frame':    365 obs. of  11 variables:
##  $ tanggal : chr  "1/1/2021" "1/2/2021" "1/3/2021" "1/4/2021" ...
##  $ pm10    : int  43 58 64 50 59 73 36 38 60 24 ...
##  $ pm25    : num  94.7 94.7 94.7 94.7 94.7 ...
##  $ so2     : int  58 86 93 67 89 81 52 68 77 39 ...
##  $ co      : int  29 38 25 24 24 29 22 26 34 16 ...
##  $ o3      : int  35 64 62 31 35 66 55 51 42 38 ...
##  $ no2     : int  65 80 86 77 77 85 72 71 80 59 ...
##  $ max     : int  65 86 93 77 89 85 72 71 80 59 ...
##  $ critical: chr  "O3" "PM25" "PM25" "O3" ...
##  $ categori: Factor w/ 3 levels "SEDANG","TIDAK SEHAT",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ location: chr  "DKI2" "DKI3" "DKI3" "DKI2" ...

Variabel categori diubah menjadi faktor dengan label baru: “SEDANG”, “TIDAK SEHAT”, dan “BAIK”. Transformasi ini penting untuk memastikan bahwa model menginterpretasikan variabel sebagai kategori nominal.

3. Multinomial Logistic Regression dengan Satu Variabel Independen (IV)

library(nnet)

# Model dengan satu IV (pm10)
fit_basic <- multinom(categori ~ pm10, data = df_MAR)

## # weights:  9 (4 variable)
## initial  value 400.993485 
## iter  10 value 141.889810
## iter  20 value 128.822646
## final  value 128.609162 
## converged

# Summary model
summary(fit_basic)

## Call:
## multinom(formula = categori ~ pm10, data = df_MAR)
## 
## Coefficients:
##             (Intercept)      pm10
## TIDAK SEHAT   -6.461683 0.2966904
## BAIK         -24.916964 0.5799524
## 
## Std. Errors:
##             (Intercept)      pm10
## TIDAK SEHAT    3.765767 0.1302711
## BAIK           4.326145 0.1344042
## 
## Residual Deviance: 257.2183 
## AIC: 265.2183

Model pertama dibuat dengan hanya menggunakan pm10 sebagai prediktor untuk memodelkan peluang relatif suatu observasi termasuk dalam masing-masing kategori kualitas udara. Hasil summary menunjukkan koefisien logit untuk masing-masing kategori dibandingkan dengan referensi dasar.

4. Output Lebih Rapi menggunakan broom dan kableExtra

library(broom)
library(kableExtra)

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

## 
## Attaching package: 'kableExtra'

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

# Tampilkan hasil model dalam bentuk tabel yang lebih rapi
tidy(fit_basic, conf.int = TRUE) %>%
  kable() %>%
  kable_styling("basic", full_width = FALSE)

y.level	term	estimate	std.error	statistic	p.value	conf.low	conf.high
TIDAK SEHAT	(Intercept)	-6.4616826	3.7657669	-1.715901	0.0861802	-13.8424500	0.9190848
TIDAK SEHAT	pm10	0.2966904	0.1302711	2.277485	0.0227573	0.0413638	0.5520170
BAIK	(Intercept)	-24.9169637	4.3261448	-5.759623	0.0000000	-33.3960516	-16.4378758
BAIK	pm10	0.5799524	0.1344042	4.314987	0.0000160	0.3165250	0.8433798

Untuk meningkatkan keterbacaan hasil model, fungsi tidy() dari paket broom digunakan. Output kemudian diformat ke dalam tabel yang lebih menarik menggunakan kableExtra.

5. Multinomial Logistic Regression dengan Multiple Variabel Independen

# Model dengan semua IV
fit_full <- multinom(categori ~ pm10 + pm25 + so2 + co + o3 + no2 + max, data = df_MAR)

## # weights:  27 (16 variable)
## initial  value 400.993485 
## iter  10 value 167.642690
## iter  20 value 53.106119
## iter  30 value 6.100912
## iter  40 value 0.518726
## iter  50 value 0.235631
## iter  60 value 0.234124
## iter  70 value 0.230342
## iter  80 value 0.229778
## iter  90 value 0.224620
## iter 100 value 0.207323
## final  value 0.207323 
## stopped after 100 iterations

# Summary model
summary(fit_full)

## Call:
## multinom(formula = categori ~ pm10 + pm25 + so2 + co + o3 + no2 + 
##     max, data = df_MAR)
## 
## Coefficients:
##              (Intercept)       pm10     pm25       so2       co        o3
## TIDAK SEHAT    0.3921943 -0.4115771 3.985658 -4.701128 1.087775 -1.438958
## BAIK        -765.0423351 -0.5473424 2.072766 -4.841885 0.882786 -1.419694
##                  no2       max
## TIDAK SEHAT 2.139859  1.037085
## BAIK        2.155423 10.730381
## 
## Std. Errors:
##             (Intercept)     pm10     pm25      so2       co       o3      no2
## TIDAK SEHAT   0.9655228 14.81667 5.996126 14.23308 2.184300 22.82341 1.405316
## BAIK          0.1045791 14.86445 9.037850 14.23878 2.196533 22.82482 1.414288
##                  max
## TIDAK SEHAT 23.36567
## BAIK         9.03785
## 
## Residual Deviance: 0.4146465 
## AIC: 32.41465

Selanjutnya, model diperluas dengan memasukkan semua polutan (pm10, pm25, so2, co, o3, no2) serta variabel max. Ini bertujuan untuk melihat bagaimana setiap variabel polutan mempengaruhi peluang relatif suatu kategori kualitas udara.

6. Output Model Full dalam Bentuk Tabel

# Menampilkan hasil model dengan semua variabel
tidy(fit_full, conf.int = TRUE) %>%
  kable() %>%
  kable_styling("basic", full_width = FALSE)

y.level	term	estimate	std.error	statistic	p.value	conf.low	conf.high
TIDAK SEHAT	(Intercept)	0.3921943	0.9655228	0.4061989	0.6845965	-1.5001956	2.284584
TIDAK SEHAT	pm10	-0.4115771	14.8166653	-0.0277780	0.9778392	-29.4517074	28.628553
TIDAK SEHAT	pm25	3.9856583	5.9961256	0.6647056	0.5062388	-7.7665319	15.737849
TIDAK SEHAT	so2	-4.7011276	14.2330751	-0.3302960	0.7411763	-32.5974421	23.195187
TIDAK SEHAT	co	1.0877749	2.1842995	0.4979972	0.6184860	-3.1933734	5.368923
TIDAK SEHAT	o3	-1.4389583	22.8234120	-0.0630475	0.9497287	-46.1720239	43.294107
TIDAK SEHAT	no2	2.1398588	1.4053160	1.5226887	0.1278366	-0.6145099	4.894228
TIDAK SEHAT	max	1.0370847	23.3656693	0.0443850	0.9645975	-44.7587856	46.832955
BAIK	(Intercept)	-765.0423351	0.1045791	-7315.4410849	0.0000000	-765.2473064	-764.837364
BAIK	pm10	-0.5473424	14.8644486	-0.0368222	0.9706267	-29.6811263	28.586441
BAIK	pm25	2.0727663	9.0378496	0.2293429	0.8186024	-15.6410935	19.786626
BAIK	so2	-4.8418848	14.2387837	-0.3400490	0.7338196	-32.7493879	23.065618
BAIK	co	0.8827860	2.1965334	0.4018997	0.6877579	-3.4223402	5.187912
BAIK	o3	-1.4196938	22.8248157	-0.0621996	0.9504039	-46.1555105	43.316123
BAIK	no2	2.1554225	1.4142879	1.5240338	0.1275003	-0.6165308	4.927376
BAIK	max	10.7303808	9.0378496	1.1872714	0.2351206	-6.9834790	28.444241

Hasil dari model multinomial penuh ditampilkan kembali dalam format tabel agar memudahkan pembacaan dan interpretasi masing-masing koefisien model.

7. Interpretasi dalam Bentuk Relative Risk Ratio (RRR)

# Mengubah koefisien logit menjadi Relative Risk Ratios
exp(coef(fit_full))

##             (Intercept)      pm10      pm25         so2       co        o3
## TIDAK SEHAT    1.480225 0.6626044 53.820710 0.009085027 2.967664 0.2371747
## BAIK           0.000000 0.5784852  7.946776 0.007892165 2.417626 0.2417880
##                  no2          max
## TIDAK SEHAT 8.498238     2.820981
## BAIK        8.631536 45724.100482

Menghitung eksponensial dari koefisien logit menghasilkan Relative Risk Ratios (RRR), yang menggambarkan perubahan faktor risiko relatif untuk masing-masing kategori dibandingkan kategori referensi.

# Menggunakan tidy() dengan exponentiate = TRUE
tidy(fit_full, conf.int = TRUE, exponentiate = TRUE) %>%
  kable() %>%
  kable_styling("basic", full_width = FALSE)

y.level	term	estimate	std.error	statistic	p.value	conf.low	conf.high
TIDAK SEHAT	(Intercept)	1.480225e+00	0.9655228	0.4061989	0.6845965	0.2230865	9.821601e+00
TIDAK SEHAT	pm10	6.626044e-01	14.8166653	-0.0277780	0.9778392	0.0000000	2.711582e+12
TIDAK SEHAT	pm25	5.382071e+01	5.9961256	0.6647056	0.5062388	0.0004237	6.836925e+06
TIDAK SEHAT	so2	9.085000e-03	14.2330751	-0.3302960	0.7411763	0.0000000	1.184518e+10
TIDAK SEHAT	co	2.967664e+00	2.1842995	0.4979972	0.6184860	0.0410332	2.146317e+02
TIDAK SEHAT	o3	2.371747e-01	22.8234120	-0.0630475	0.9497287	0.0000000	6.344420e+18
TIDAK SEHAT	no2	8.498238e+00	1.4053160	1.5226887	0.1278366	0.5409059	1.335168e+02
TIDAK SEHAT	max	2.820981e+00	23.3656693	0.0443850	0.9645975	0.0000000	2.184208e+20
BAIK	(Intercept)	0.000000e+00	0.1045791	-7315.4410849	0.0000000	0.0000000	0.000000e+00
BAIK	pm10	5.784852e-01	14.8644486	-0.0368222	0.9706267	0.0000000	2.599763e+12
BAIK	pm25	7.946776e+00	9.0378496	0.2293429	0.8186024	0.0000002	3.919426e+08
BAIK	so2	7.892200e-03	14.2387837	-0.3400490	0.7338196	0.0000000	1.040569e+10
BAIK	co	2.417626e+00	2.1965334	0.4018997	0.6877579	0.0326360	1.790943e+02
BAIK	o3	2.417880e-01	22.8248157	-0.0621996	0.9504039	0.0000000	6.485644e+18
BAIK	no2	8.631536e+00	1.4142879	1.5240338	0.1275003	0.5398139	1.380169e+02
BAIK	max	4.572410e+04	9.0378496	1.1872714	0.2351206	0.0009271	2.255157e+12

Dalam output ini, RRR dan interval kepercayaannya (confidence interval) juga ditampilkan agar lebih mudah diinterpretasikan. Nilai RRR > 1 menunjukkan peningkatan risiko relatif, sementara RRR < 1 menunjukkan penurunan.

8. Marginal Effects Analysis

library(marginaleffects)

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

# Marginal effects untuk pm10
mfx_pm10 <- avg_comparisons(fit_full, variables = "pm10", type = "probs")
mfx_pm10

## 
##        Group  Estimate Std. Error       z Pr(>|z|)   S     2.5 %   97.5 %
##  SEDANG       3.06e-06   0.000129  0.0236    0.981 0.0 -0.000251 0.000257
##  TIDAK SEHAT  7.08e-05   0.000680  0.1041    0.917 0.1 -0.001261 0.001403
##  BAIK        -7.38e-05   0.000667 -0.1107    0.912 0.1 -0.001381 0.001234
## 
## Term: pm10
## Type:  probs 
## Comparison: +1

# Marginal effects untuk co
mfx_co <- avg_comparisons(fit_full, variables = "co", type = "probs")
mfx_co

## 
##        Group  Estimate Std. Error       z Pr(>|z|)   S     2.5 %   97.5 %
##  SEDANG      -1.01e-05   0.000209 -0.0484    0.961 0.1 -0.000419 0.000399
##  TIDAK SEHAT  1.22e-04   0.000711  0.1715    0.864 0.2 -0.001272 0.001515
##  BAIK        -1.12e-04   0.000672 -0.1664    0.868 0.2 -0.001429 0.001205
## 
## Term: co
## Type:  probs 
## Comparison: +1

Analisis efek marginal dilakukan untuk mengukur seberapa besar perubahan probabilitas kategori akibat perubahan kecil pada variabel prediktor (pm10 dan co). Ini memberikan interpretasi yang lebih intuitif tentang hubungan antara prediktor dan outcome dibandingkan interpretasi berbasis log-odds.

Kesimpulan

Melalui analisis multinomial logistic regression ini, diperoleh gambaran bagaimana kadar polutan seperti pm10, pm25, so2, co, o3, no2, dan max mempengaruhi kualitas udara yang dikategorikan sebagai “SEDANG”, “TIDAK SEHAT”, dan “BAIK”. Koefisien logit dan relative risk ratios memberikan informasi arah dan kekuatan hubungan. Sementara itu, analisis efek marginal membantu menginterpretasikan pengaruh masing-masing variabel dalam satuan probabilitas, yang lebih mudah dipahami oleh pengambil keputusan.