MODEL REGRESI LOGISTIK MULTINOMIAL

Tujuan

  1. Mempelajari program R untuk membuat model regresi logistic multinomial.
  2. Mengaplikasikan program R untuk membuat model regresi logistic multinomial.

Pendahuluan

Regresi logistik multinomial adalah model regresi yang digunakan untuk menyelesaikan kasus regresi dengan variabel dependen berupa data kualitatif berbentuk multinomial (lebih dari dua kategori) dengan satu atau lebih variabel independen. Persamaan model regresi logistic multinomial dapat dituliskan sebagai berikut: \(g_j\) \((x)=β_{j0}+β_{j1} x_1+β_{j2} x_{2}+⋯+β_{jp} x_{p}\), dengan \(g_j (x)\) merupakan variable respon yang berupa variable kategori dengan skala pengukuran nominal, \(x_{p}\) menyatakan variable independent, dan \(β_{jp}\) adalah parameter.

Rumus dasar model persamaan regresi logistik multinomial

Metode yang digunakan untuk mengestimasi parameter model regresi logistik multinomial pada penulisan ini adalah metode maksimum likelihood (maximum likelihood methods). Persamaan likelihood pada regresi logistik multinomial merupakan persamaan nonlinear dalam parameter koefisien regresi \(β_jp\), sehingga untuk menyelesaikan persamaan tersebut sampai diperoleh nilai estimasi parameternya digunakan algoritma Newton Raphson. Kemudian setelah diperoleh estimasi parameter, dilakukan uji taraf nyata parameter menggunakan uji Rasio likelihood dan uji Wald

Program R

Contoh Kasus :
Kumpulan data berisi variabel pada 200 siswa. Dimana Prog (program type) adalah variable dependen (Y) dan Ses (Social Economic status) dan Write (Writing Score) adalah variable independen (\(X_1\) dan \(X_2\))

Pagkage

library(foreign)
library(nnet)
library(ggplot2)
library(reshape2)

Pemanggilan Data

ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")

Membuat tabel kontingensi dari Y dan X

with(ml, table(ses, prog))
##         prog
## ses      general academic vocation
##   low         16       19       12
##   middle      20       44       31
##   high         9       42        7

Estimasi model regresi logistik multinomial

ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)
## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.982880
## final  value 179.981726 
## converged
summary(test)
## Call:
## multinom(formula = prog2 ~ ses + write, data = ml)
## 
## Coefficients:
##          (Intercept)  sesmiddle    seshigh      write
## general     2.852198 -0.5332810 -1.1628226 -0.0579287
## vocation    5.218260  0.2913859 -0.9826649 -0.1136037
## 
## Std. Errors:
##          (Intercept) sesmiddle   seshigh      write
## general     1.166441 0.4437323 0.5142196 0.02141097
## vocation    1.163552 0.4763739 0.5955665 0.02221996
## 
## Residual Deviance: 359.9635 
## AIC: 375.9635

Output ringkasan model memiliki blok koefisien dan blok kesalahan standar. Masing-masing blok memiliki satu baris nilai yang sesuai dengan persamaan model. Pada blok koefisien, baris pertama membandingkan prog = “general” dengan baseline prog = “academic”, dan baris kedua membandingkan prog = “vocation” dengan baseline prog = “academic”. Dari koefisien, kita bisa menuliskan persamaan model untuk kedua kategori tersebut.

\(β_{13}\) : Peningkatan satu unit pada variabel write dikaitkan dengan penurunan log odds berada di program general dibandingkan dengan academic program sebesar 0.058.
\(β_{23}\) : Peningkatan satu unit pada variabel write dikaitkan dengan penurunan log odds berada di program vocation dibandingkan dengan academic program sebesar 0.1136.
\(β_{12}\) : Log odds berada di program general vs. academic akan menurun sebesar 1.163 jika berpindah dari SES “low” ke “high”.
\(β_{11}\) : Log odds berada di program general vs. academic akan menurun sebesar 0.533 jika berpindah dari SES “low” ke “middle”. Meskipun tidak signifikan
\(β_{22}\) : Log odds berada di program vocation vs. academic akan menurun sebesar 0.983 jika berpindah dari SES “low” ke “high”.
\(β_{21}\) : Log odds berada di program vocation vs. academic akan menurun sebesar 0.291 jika berpindah dari SES “low” ke “middle”. Meskipun tidak signifikan

Hitung nilai z untuk koefisien

z <- summary(test)$coefficients/summary(test)$standard.errors
z
##          (Intercept)  sesmiddle   seshigh     write
## general     2.445214 -1.2018081 -2.261334 -2.705562
## vocation    4.484769  0.6116747 -1.649967 -5.112689

Uji wald untuk menghitung p-value

p <- (1 - pnorm(abs(z), 0, 1)) * 2
p
##           (Intercept) sesmiddle    seshigh        write
## general  0.0144766100 0.2294379 0.02373856 6.818902e-03
## vocation 0.0000072993 0.5407530 0.09894976 3.176045e-07

Odds Ration

Esponensialkan koefisien untuk mendapatkan odds ratio

exp(coef(test))
##          (Intercept) sesmiddle   seshigh     write
## general     17.32582 0.5866769 0.3126026 0.9437172
## vocation   184.61262 1.3382809 0.3743123 0.8926116

Ekspansi koefisien log-odds ke dalam bentuk odds ratio memungkinkan kita untuk menginterpretasikan dampak dari variabel prediktor pada probabilitas relatif. Ini membantu memahami seberapa besar peluang relatif perubahan hasil ketika variabel prediktor berubah, dengan semua variabel lainnya tetap konstan. Yang dimana berdasarkan hasil sintak
• Rasio risiko relatif untuk peningkatan satu unit pada variabel write adalah 0.9437 untuk berada di program general dibandingkan dengan academic.
• Rasio risiko relatif berpindah dari ses = 1 ke ses = 3 adalah 0.3126 untuk berada di program general dibandingkan academic.

Prediksi probabilitas menggunakan model

pp <- fitted(test)
head (pp)
##    academic   general  vocation
## 1 0.1482764 0.3382454 0.5134781
## 2 0.1202017 0.1806283 0.6991700
## 3 0.4186747 0.2368082 0.3445171
## 4 0.1726885 0.3508384 0.4764731
## 5 0.1001231 0.1689374 0.7309395
## 6 0.3533566 0.2377976 0.4088458

Prediksi Porbabilitas untuk data baru

dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
predict(test, newdata = dses, "probs")
##    academic   general  vocation
## 1 0.4396845 0.3581917 0.2021238
## 2 0.4777488 0.2283353 0.2939159
## 3 0.7009007 0.1784939 0.1206054

Fungsi ini digunakan untuk menghitung prediksi probabilitas dari model regresi logistik multinomial untuk setiap kategori hasil. Dengan menggunakan fungsi fitted(), kamu dapat menghasilkan probabilitas prediksi untuk setiap observasi di dataset. Ini membantu memahami bagaimana model memprediksi hasil berdasarkan variabel input. Selain itu, dengan menggunakan dataset kecil seperti dses, kamu bisa melihat bagaimana perubahan variabel (misalnya ses) memengaruhi probabilitas prediksi dengan variabel lain (seperti write) yang tetap konstan. Ini membantu analisis sensitivitas model terhadap variabel prediktor.

Melihat probabilitas Prediksi rata-rata
Untuk berbagai nilai variabel prediktor kontinu write di dalam setiap level ses (status sosial ekonomi)

dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),
                                                                                   3))
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))
by(pp.write[, 3:5], pp.write$ses, colMeans)
## pp.write$ses: high
##  academic   general  vocation 
## 0.6164315 0.1808037 0.2027648 
## ------------------------------------------------------------ 
## pp.write$ses: low
##  academic   general  vocation 
## 0.3972977 0.3278174 0.2748849 
## ------------------------------------------------------------ 
## pp.write$ses: middle
##  academic   general  vocation 
## 0.4256198 0.2010864 0.3732938

Maksud dari sintaks ini adalah untuk memahami model regresi logistik multinomial dengan melihat probabilitas prediksi rata-rata untuk berbagai nilai variabel prediktor kontinu write di dalam setiap level ses (status sosial ekonomi). Ini dilakukan dengan cara: 1. Membuat dataset untuk setiap kombinasi nilai ses dan write. 2. Menghitung probabilitas prediksi untuk setiap kombinasi tersebut. 3. Mengambil rata-rata probabilitas dalam setiap level ses. Hasil ini membantu dalam menganalisis bagaimana nilai write memengaruhi probabilitas hasil pada setiap level ses.

Plot
Untuk melihat grafiknya

lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")
head(lpp)  # view first few rows
##   ses write variable probability
## 1 low    30 academic  0.09843588
## 2 low    31 academic  0.10716868
## 3 low    32 academic  0.11650390
## 4 low    33 academic  0.12645834
## 5 low    34 academic  0.13704576
## 6 low    35 academic  0.14827643
ggplot(lpp, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~
                                                                                        ., scales = "free")

Kegunaan dari visualisasi ini adalah untuk membantu memahami bagaimana probabilitas prediksi berubah terhadap nilai variabel write pada setiap level ses (status sosial ekonomi) untuk setiap kategori hasil (akademik, umum, panggilan vokasi). Plot ini menunjukkan bagaimana skor menulis memengaruhi peluang berada dalam setiap program untuk status sosial ekonomi yang berbeda. Visualisasi ini membantu menganalisis pola prediksi yang dihasilkan oleh model regresi logistik multinomial secara lebih intuitif dan mendalam.

Untuk Latihan silahkan gunakan yang ada pada modul kawan-kawan yaaa