Analisis Regresi Logistik

Alfidhia

2022-11-19

library(car)

Data

Data yang digunakan adalah data Kecelakaan Kapal Titanic pada tanggal 15 April 1912. Data ini berisi 9 variabel dengan 889 observasi/amatan. Namun, banyak variabel yang digunakan untuk analisis ini hanyalah 4, yaitu:

  1. Survived : Selamat/tidaknya penumpang (Biner)

  2. Pclass : Kelas masing-masing observasi/penumpang (Kategorik/Ordinal)

  3. Sex : Jenis kelamin penumpang (Biner)

  4. Age : Umur penumpang (Numerik/Rasio)

Secara default, R akan membaca peubah Survived dan Pclass sebagai data numerik dan peubah Sex sebagai peubah karakter. Syntax dibawah dijalankan untuk menyesuaikan jenis data yang seharusnya.

titanic <- read.csv("D:/Semester 5/4.1. Pengantar Analisis Data Kategorik/titanic_clean.csv")
titanic <- titanic[2:5]
titanic$Survived <- as.factor(titanic$Survived)
titanic$Pclass <- as.factor(titanic$Pclass)
titanic$Sex <- as.factor(titanic$Sex)
str(titanic)
## 'data.frame':    889 obs. of  4 variables:
##  $ Survived: Factor w/ 2 levels "0","1": 1 2 2 2 1 1 1 1 2 2 ...
##  $ Pclass  : Factor w/ 3 levels "1","2","3": 3 1 3 1 3 3 1 3 3 2 ...
##  $ Sex     : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...
##  $ Age     : num  22 38 26 35 35 ...

Dari output di atas, dapat terlihat bahwa secara default R memlih referensi peubah Survived adalah 0 (tidak selamat), referensi peubah Pclass adalah 1, dan referensi peubah Sex adalah “female”.

Analisis Regresi Logistik

Model Regresi Logistik

Model pada kasus ini adalah sebagai berikut:

\[ logit[\pi(x)]=\alpha+\beta_1c_1+\beta_2c_2+\beta_3d+\beta_4x \]

dengan \(\alpha\) sebagai intersep, \(\beta\) sebagai koefisien tiap peubah, \(c\) sebagai peubah Pclass, \(c_1\) bernilai 1 untuk Pclass2, bernilai 0 lainnya dan \(c_2\) bernilai 1 untuk Pclass3, bernilai 0 lainnya, \(d\) sebagai peubah Sex, bernilai 1 untuk laki-laki, bernilai 0 lainnya, dan \(x\) sebagai peubah Age.

Pendugaan Parameter

reg_log <- glm(Survived ~ ., data = titanic, family = binomial(link = "logit"))
summ <- summary(reg_log)
summ
## 
## Call:
## glm(formula = Survived ~ ., family = binomial(link = "logit"), 
##     data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6469  -0.6630  -0.4203   0.6292   2.4285  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.539586   0.365376   9.688  < 2e-16 ***
## Pclass2     -1.115005   0.257824  -4.325 1.53e-05 ***
## Pclass3     -2.321885   0.240989  -9.635  < 2e-16 ***
## Sexmale     -2.603906   0.186871 -13.934  < 2e-16 ***
## Age         -0.033530   0.007385  -4.540 5.62e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1182.82  on 888  degrees of freedom
## Residual deviance:  804.68  on 884  degrees of freedom
## AIC: 814.68
## 
## Number of Fisher Scoring iterations: 5

Dari output di atas, dapat terlihat bahwa setiap peubah berpengaruh terhadap model. Sehingga didapat modelnya yaitu:

\[ logit[\hat{\pi(x)}]=3.539-1.115c_1-2.321c_2-2.6d-0.033x \]

Pendugaan peluang peubah kategorik

Untuk menduga peluang peubah kategorik, peubah numerik akan dianggap konstan. Dalam hal ini, kita akan menganggap peubah Age konstan di nilai rata-ratanya, yaitu sebesar 29.65345 tahun. Nilai ini didapatkan dengan syntax:

avg_age <- mean(titanic$Age)
avg_age
## [1] 29.65345

Peluang peubah kategorik didapatkan dengan fungsi berikut:

\[ \hat{\pi(x)}=\frac{exp(logit[\hat{\pi(x)}]}{1+exp(logit[\hat{\pi(x)}])} \]

dapat dibuat fungsi di r sebagai berikut:

p <- function(logit){
  exp(logit)/(1+exp(logit))
}

dengan fungsi logitnya:

alp <- summ$coefficients[1,1]
b1 <- summ$coefficients[2,1]
b2 <- summ$coefficients[3,1]
b3 <- summ$coefficients[4,1]
b4 <- summ$coefficients[5,1]
logit_x <- function(c1,c2,d){
  alp+b1*c1+b2*c2+b3*d+b4*avg_age
}

sehingga, dapat didapatkan peluang masing-masing peubah kategoriknya:

peubah=c("Pclass1_Pria","Pclass1_Wanita","Pclass2_Pria","Pclass2_Wanita","Pclass3_Pria","Pclass3_Wanita")
logit=c(logit_x(0,0,1),logit_x(0,0,0),logit_x(1,0,1),logit_x(1,0,0),logit_x(0,1,1),logit_x(0,1,0))
peluang=c(p(logit_x(0,0,1)),p(logit_x(0,0,0)),p(logit_x(1,0,1)),p(logit_x(1,0,0)),p(logit_x(0,1,1)),p(logit_x(0,1,0)))
keselamatan=ifelse(peluang>0.5,"Selamat","Tidak selamat")
data.frame(peubah, logit, peluang, keselamatan)

Dari tabel di atas, dapat terlihat bahwa pada umur 29 tahun-an, perempuan cenderung berpeluang untuk selamat jika dibanding laki-laki. Selain itu, dapat terlihat juga bahwa semakin besar kelas P, maka peluang keselamatan penumpang semakin rendah.

Perhitungan Odds Rasio

Data Kategorik

Untuk dapat mengetahui perbandingan antara peubah kategorik lain dengan peubah kategorik referensi, kita perlu menghitung nilai eksponensial dari masing-masing parameter peubah dummy-nya. Misalnya untuk membandingkan dugaan keselamatan male dibanding female, kita dapat mengeksponensialkan nilai \(\beta_3\) yang merupakan koefisien bagi peubah dummy Sex. Cara tadi diulang pada setiap peubah dummy sehingga didapatkan dua tabel sebagai berikut:

Male vs Female

data.frame(Gender="Laki-laki",Odds_Rasio=exp(b3))

Interpretasi: Untuk nilai peubah lainnya tetap, dugaan odds bagi seorang laki-laki untuk selamat adalah 0.07 kali dibanding perempuan.

Pclass: 1 vs 2 dan 1 vs 3

data.frame(Pclass=c("Pclass2", "Pclass3"),Odds_Rasio=c(exp(b1),exp(b2)))

  • Interpretasi Pclass 2: Untuk nilai peubah lainnya tetap, dugaan odds bagi orang yang berada di kelas 2 untuk selamat adalah 0.32 kali dibanding orang yang berada di kelas 1.

  • Interpretasi Pclass 3: Untuk nilai peubah lainnya tetap, dugaan odds bagi orang yang berada di kelas 3 untuk selamat adalah 0.09 kali dibanding orang yang berada di kelas 1.

Data Numerik

Age: x vs x+1

Untuk dapat membandingkan nilai odds rasio data numerik, dalam hal ini Age, kita dapat membandingkan dua nilai Age dengan beda satu. Fungsi yang digunakan adalah sebagai berikut:

\[ exp(\beta_4)=OR=\frac{Odds(X=x+1)}{Odds(X=x)} \]

data.frame("Age"=c("x+1"), Odds_Rasio=exp(b4))

Interpretasi: Peluang keselamatan penumpang akan menurun secara multiplikatif sebesar 0.033 atau 3.3% untuk setiap kenaikan 1 tahun usia penumpang.

Uji Parsial - Wald

\[ \text{Hipotesis:}\\ H_0:\beta_j=0\\ H_1:\beta_j\neq0 \]

Anova(reg_log, type="II", test="Wald")

Berdasarkan hasil uji di atas, dapat disimpulkan bahwa setiap peubah yaitu Pclass, sex, dan age berpengaruh signifikan terhadap keselematan penumpang kapal titanic pada taraf nyata 0.1%.

Uji Simultan - Likelihood Rasio

\[ \text{Hipotesis:}\\ H_0:\beta_1=\beta_2=\beta_3=\beta_4=0\\ H_1:\text{minimal ada satu }\beta_j\neq0 \]

anova(reg_log,update(reg_log, ~1),test="Chisq")

Berdasarkan hasil uji di atas, dapat disimpulkan bahwa secara bersama-sama, peubah Pclass, sex, dan age berpengaruh signifikan terhadap keselematan penumpang kapal titanic pada taraf nyata 0.1%.