Project ini merupakan Project Kelompok Mata Kuliah Analisis Data Kategorik
Kapal Titanic merupakan kapal laut terbesar dan termewah yang pernah dibangun. Panjangnya 833 kaki atau sekitar 254 meter dari buritan ke busur. Lambungnya dibagi menjadi 16 kompartemen dan kapal ini dianggap kedap air. Dianggap kedap air karena sebanyak 4 kompartemen dapat dibanjiri tanpa menyebabkan hilangnya daya apung kapal. Itu juga membuat orang mengira Titanic tidak dapat tenggelam.
Namun siapa sangka tepat sebelum tengah malam pada 14 April, Titanic gagal mengalihkan jalurnya dari gunung es dan menghancurkan setidaknya 5 kompartemen lambungnya. Karena kompartemen Titanic tidak ditutup di bagian atas, air dari kompartemen yang pecah memenuhi setiap kompartemen. Itu menyebabkan busur tenggelam dan buritan naik hingga hampir tegak lurus. Setelah itu Titanic pecah menjadi 2 bagian. Begitulah hingga kapal tenggelam ke dasar samudera. Maka dari itu kami tertarik untuk membahas lebih lanjut faktor-faktor apa saja yang mempengaruhi sehingga penumpang dapat bertahan hidup (survive) menggunakan analisis regresi logistik. Regresi logistik adalah sebuah pendekatan untuk membuat model prediksi seperti halnya regresi linear atau yang biasa disebut dengan istilah Ordinary Least Squares (OLS) regression.
Perbedaannya adalah pada regresi logistik, peneliti memprediksi variabel terikat yang berskala dikotomi. Skala dikotomi yang dimaksud adalah skala data nominal dengan dua kategori. Karena dari data yang kami peroleh, kebanyakan variabelnya berskala dikotomi, mulai dari survival, kelas, jenis kelamin, dan titik mulai (embarked).
Data yang digunakan dalam penelitian ini adalah data data titanic dalam file format csv digunakan untuk menganalisis faktor-faktor yang mempengaruhi penumpang bertahan hidup atau tidak (survived) dengan jumlah observasi yaitu terhadap 712 orang. Secara total data terdiri dari 11 variabel:
Pengolahan data penelitian ini dibantu dengan menggunakan Software Rstudio. Langkah-langkah pengolahan data adalah sebagai berikut:
Data Titanic diinput untuk dilakukan penelitian dengan 11 variabel. Akan digunakan 1 variabel sebagai variabel dependent dan variabel lainnya akan diindikasi sebagai variabel yang mungkin menjadi factor atau variabel independent.
MyData <- read.csv(file.choose(), header=TRUE)
head(MyData)
## PassengerId Survived Pclass
## 1 1 0 3
## 2 2 1 1
## 3 3 1 3
## 4 4 1 1
## 5 5 0 3
## 6 7 0 1
## Name Sex Age SibSp Parch
## 1 Braund, Mr. Owen Harris male 22 1 0
## 2 Cumings, Mrs. John Bradley (Florence Briggs Thayer) female 38 1 0
## 3 Heikkinen, Miss. Laina female 26 0 0
## 4 Futrelle, Mrs. Jacques Heath (Lily May Peel) female 35 1 0
## 5 Allen, Mr. William Henry male 35 0 0
## 6 McCarthy, Mr. Timothy J male 54 0 0
## Ticket Fare Cabin Embarked
## 1 A/5 21171 7.2500 S
## 2 PC 17599 71.2833 C85 C
## 3 STON/O2. 3101282 7.9250 S
## 4 113803 53.1000 C123 S
## 5 373450 8.0500 S
## 6 17463 51.8625 E46 S
Survived <- MyData$Survived
levels (Survived) = c("Not Survived","Survived")
Sex <- as.factor(MyData$Sex)
Age <- MyData$Age
MyData$Pclass <- MyData$Pclass-1
MyData$Pclass <- as.factor(MyData$Pclass)
Fare <- MyData$Fare
Ticket <- MyData$Ticket
Sibling <- MyData$SibSp
Parents <- MyData$Parch
MyData$Embarked <- as.factor(MyData$Embarked)
Pemilihan Model dengan menguji variabel independent dan mencari 3 variabel independent yang paling berpengaruh dilakukan dengan menggunakan variabel Sex, Age, Passanger Class (Pclass), Passanger Fare (Fare), Number of Siblings/Sibsp (Sibling), Number of Parents/parch (Parents) yang dimungkinkan menjadi faktor dari variabel dependent Survived dengan Sex, Passanger Class, Embarked, dan Survived dijadikan sebagai variabel kategorik. Dengan bantuan software R diperoleh nilai sebagai berikut :
Titanic.logit <- glm(Survived ~ Sex + Age + Pclass + Fare + Sibling + Parents , data = MyData , family = binomial)
summary(Titanic.logit)
##
## Call:
## glm(formula = Survived ~ Sex + Age + Pclass + Fare + Sibling +
## Parents, family = binomial, data = MyData)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7984 -0.6437 -0.3886 0.6353 2.4422
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.172041 0.503349 8.289 < 2e-16 ***
## Sexmale -2.628356 0.220355 -11.928 < 2e-16 ***
## Age -0.044334 0.008276 -5.357 8.46e-08 ***
## Pclass1 -1.285746 0.321735 -3.996 6.43e-05 ***
## Pclass2 -2.495085 0.338640 -7.368 1.73e-13 ***
## Fare 0.002022 0.002557 0.791 0.42904
## Sibling -0.375639 0.127328 -2.950 0.00318 **
## Parents -0.059443 0.122898 -0.484 0.62861
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 960.90 on 711 degrees of freedom
## Residual deviance: 635.27 on 704 degrees of freedom
## AIC: 651.27
##
## Number of Fisher Scoring iterations: 5
Residual Deviance = 635.27 adalah deviance dari model dengan variabel independent Sex, Age, Pclass, Fare, Sibling, Parents sebagai variabel independent. 𝐻1 = logit (𝑃(𝑌 = 1)) = 𝛼 + 𝛽1 Sex + 𝛽2 Age + 𝛽3 Pclass + 𝛽4 Fare + 𝛽5 Sibling + 𝛽6 Parents
Null Deviance = 960.90 adalah deviance dari model tanpa variabel independent. 𝐻0 = logit (𝑃(𝑌 = 1)) = α
Dapat dilihat output diatas dari p-value yang lebih kecil dari alpha (0.05) adalah variable yang signifikan yaitu variable Sexmale, Age, Pclass2, Pclass3, dan Sibling. Dalam penelitian ini kita mengkhususkan 3 variabel independent yang berpengaruh terhadap survived dan dipilih variabel yang menunujukan paling signifikan yaitu variabel Sex, Age dan Pclass.
Selanjutnya melakukan seleksi model dengan eliminasi backward dimulai dengan model yang kompleks ke model yang lebih sederhana, dengan membuang variabel yang paling tidak signifikan.
Pertama-tama cek siginifikansi interaksi 3 arah, dengan model yang melihabatkan interaksi 3 arah dan di dapatkan nilai P-value sebagai berikut.
Titanic.logitmodel <- glm(Survived ~ Age*Sex*Pclass , data = MyData , family = binomial)
step(Titanic.logitmodel, test ="Chisq")
## Start: AIC=625.27
## Survived ~ Age * Sex * Pclass
##
## Df Deviance AIC LRT Pr(>Chi)
## - Age:Sex:Pclass 2 604.18 624.18 2.9103 0.2334
## <none> 601.27 625.27
##
## Step: AIC=624.18
## Survived ~ Age + Sex + Pclass + Age:Sex + Age:Pclass + Sex:Pclass
##
## Df Deviance AIC LRT Pr(>Chi)
## <none> 604.18 624.18
## - Age:Sex 1 606.70 624.70 2.5148 0.11278
## - Age:Pclass 2 610.90 626.90 6.7167 0.03479 *
## - Sex:Pclass 2 631.67 647.67 27.4862 1.075e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call: glm(formula = Survived ~ Age + Sex + Pclass + Age:Sex + Age:Pclass +
## Sex:Pclass, family = binomial, data = MyData)
##
## Coefficients:
## (Intercept) Age Sexmale Pclass1
## 3.403838 -0.003488 -2.469639 1.260834
## Pclass2 Age:Sexmale Age:Pclass1 Age:Pclass2
## -3.162357 -0.029972 -0.064136 -0.014957
## Sexmale:Pclass1 Sexmale:Pclass2
## -1.484805 1.667659
##
## Degrees of Freedom: 711 Total (i.e. Null); 702 Residual
## Null Deviance: 960.9
## Residual Deviance: 604.2 AIC: 624.2
Dapat dilihat dari p-value yang didapat = 0.2334 yang artinya lebih besar dari alpha maka interaksi 3 arah tidak signifikan dan selanjutnya interaksi tersebut dibuang. Setelah membuang interaksi 3 arah, sekarang kita lihat pula jika model hanya melibatkan interaksi 2 arah.
Titaniclogit <- glm(Survived ~ Age*Sex + Age*Pclass + Sex*Pclass , data = MyData , family = binomial)
drop1(Titaniclogit , test = "Chisq")
## Single term deletions
##
## Model:
## Survived ~ Age * Sex + Age * Pclass + Sex * Pclass
## Df Deviance AIC LRT Pr(>Chi)
## <none> 604.18 624.18
## Age:Sex 1 606.70 624.70 2.5148 0.11278
## Age:Pclass 2 610.90 626.90 6.7167 0.03479 *
## Sex:Pclass 2 631.67 647.67 27.4862 1.075e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Diantara semua interaksi 2 arah, didapat Age * Sex yang memiliki p-value terbesar atau Age * Sex paling tidak signifikan, sehingga Age*Sex dieliminasi.
Setelah membuang interaksi 2 arah dari Age * Sex, sekarang kita membuat model yang hanya melihabatkan interaksi 2 arah yang melibatkan Age * Pclass dan Pclass*Sex dan di dapatkan p-value sebagai berikut:
fit1 <- glm(Survived ~ Age*Pclass + Sex*Pclass , data = MyData , family = binomial)
drop1(fit1 , test = "Chisq")
## Single term deletions
##
## Model:
## Survived ~ Age * Pclass + Sex * Pclass
## Df Deviance AIC LRT Pr(>Chi)
## <none> 606.70 624.70
## Age:Pclass 2 613.16 627.16 6.463 0.0395 *
## Pclass:Sex 2 645.09 659.09 38.393 4.603e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Model dipilih karena pada interaksi 2 arah yang melibatkan AgePclass dan PclassSex parameter sudah signifikan dan p-value < α = 0,05. Dengan nama model fit2 sebagai berikut.
fit2 <- glm(Survived ~ Age + Pclass + Sex + Age*Pclass + Sex*Pclass , data = MyData , family = binomial)
fit2$coefficients
## (Intercept) Age Pclass1 Pclass2 Sexmale
## 4.40487096 -0.03049152 1.03846951 -3.84339055 -3.58911246
## Age:Pclass1 Age:Pclass2 Pclass1:Sexmale Pclass2:Sexmale
## -0.05764920 -0.00295660 -1.34616341 2.12648080
logit (𝜋) = 4,40487 − 0,03049152(Age) + 1,03846951(Pclass2) − 3,84339055(Pclass3) − 3,58911246(Sexmale) − 0,0576492(Age * Pclass2) − 0,0029566(Age * Pclass3) − 1,34616341(Pclass2 * Sexmale) + 2,12648080(Pclass3 * Sexmale)
Akan dilihat pula jika Age*Pclass di eliminasi dan diberi nama model fit3 sebagai berikut.
fit3 <- glm(Survived ~ Age + Pclass + Sex + Sex*Pclass , data = MyData , family = binomial)
fit3$coefficients
## (Intercept) Age Pclass1 Pclass2 Sexmale
## 4.87000167 -0.04205321 -1.11651587 -4.12683112 -3.59466605
## Pclass1:Sexmale Pclass2:Sexmale
## -0.71448300 2.14828534
logit(𝜋) = 4,87000167 − 0,04205321(Age) − 1,11651587(Pclass2) − 4,12683112(Pclass3) − 3,59466605(Sexmale) − 0,71448300(Pclass2 * Sexmale) + 2,14828534(Pclass3 * Sexmale)
Untuk membandingkan mana kan model yang lebih fit dapat dilakukan dengan berbagai cara:
Disini kita akan membandingan nilai AIC yang terdapat di mode fit2 dan model fit3.
fit2
##
## Call: glm(formula = Survived ~ Age + Pclass + Sex + Age * Pclass +
## Sex * Pclass, family = binomial, data = MyData)
##
## Coefficients:
## (Intercept) Age Pclass1 Pclass2
## 4.404871 -0.030492 1.038470 -3.843391
## Sexmale Age:Pclass1 Age:Pclass2 Pclass1:Sexmale
## -3.589112 -0.057649 -0.002957 -1.346163
## Pclass2:Sexmale
## 2.126481
##
## Degrees of Freedom: 711 Total (i.e. Null); 703 Residual
## Null Deviance: 960.9
## Residual Deviance: 606.7 AIC: 624.7
fit3
##
## Call: glm(formula = Survived ~ Age + Pclass + Sex + Sex * Pclass, family = binomial,
## data = MyData)
##
## Coefficients:
## (Intercept) Age Pclass1 Pclass2
## 4.87000 -0.04205 -1.11652 -4.12683
## Sexmale Pclass1:Sexmale Pclass2:Sexmale
## -3.59467 -0.71448 2.14829
##
## Degrees of Freedom: 711 Total (i.e. Null); 705 Residual
## Null Deviance: 960.9
## Residual Deviance: 613.2 AIC: 627.2
Didapat dilihat dari hasil output diatas nilai AIC lebih rendah/model lebih baik terdapat di model fit2.
Maka disini kita akan membandingkan nilai aktual dan nilai prediksi yang terdapat pada model fit2 dan model fit3.
#model fit2
prob.prediksi <- predict(fit2 ,type = "response")
table(MyData$Survived,prob.prediksi>0.5)
##
## FALSE TRUE
## 0 393 31
## 1 106 182
#Akurasi model fit2-> 393+182/712 = 0,8075 atau 80%
#model fit2
prob.prediksi2 <- predict(fit3 ,type = "response")
table(MyData$Survived,prob.prediksi2 >0.5)
##
## FALSE TRUE
## 0 386 38
## 1 108 180
#Akurasi model fit3-> 386+180/712 = 0,79494 atau 79%
Dapat dilihat dari perhitungan akurasi untuk model fit2 dan perhitungan akurasi untuk model fit3 terlihat bahwa lebih besar perhitungan akurasi untuk model fit2 maka didapat model terbaik yaitu model fit2.
Maka disini kita akan melakukan uji Hosmer and Lesmeshow goodness (hoslem.test) pada model fit2. Dengan, 𝐻0 = Model dapat diterima 𝐻1 = Model tidak dapat diterima Hasil yang didapat sebagai berikut.
library(ResourceSelection)
## ResourceSelection 0.3-5 2019-07-22
hoslem.test(Survived,fitted(fit2))
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: Survived, fitted(fit2)
## X-squared = 7.6642, df = 8, p-value = 0.4669
Dapat dilihat dari tabel diatas didapat p-value = 0.4669 > 𝛼 = 0.05 maka 𝐻0 diterima atau model dapat diterima. Dan didapatkan model terpilih yaitu Age + Pclass + Sex + Age * Pclass + Sex * Pclass.
logit (𝜋) = 4,40487 − 0,03049152(Age) + 1,03846951(Pclass2) − 3,84339055(Pclass3) − 3,58911246(Sexmale) − 0,0576492(Age * Pclass2) − 0,0029566(Age * Pclass3) − 1,34616341(Pclass2 * Sexmale) + 2,12648080(Pclass3 * Sexmale)
dengan akurasi model 0,8075 atau 81%
Untuk memudahkan dalam interpretasi, maka variabel Pclass dapat di Simplifikasi dengan menjadikannya variabel kuantitatif. Diperoleh model yaitu Age + Pclass + Sexmale + Age * Pclass + Pclass * Sexmale
MyData$Pclass <- as.numeric(as.character(MyData$Pclass))
fit2 <- glm(Survived ~ Age + Pclass + Sex + Age*Pclass + Sex*Pclass , data = MyData , family = binomial)
fit2
##
## Call: glm(formula = Survived ~ Age + Pclass + Sex + Age * Pclass +
## Sex * Pclass, family = binomial, data = MyData)
##
## Coefficients:
## (Intercept) Age Pclass Sexmale Age:Pclass
## 5.400099 -0.034001 -2.293156 -4.645139 -0.004016
## Pclass:Sexmale
## 1.488765
##
## Degrees of Freedom: 711 Total (i.e. Null); 706 Residual
## Null Deviance: 960.9
## Residual Deviance: 623.7 AIC: 635.7
prob.prediksi <- predict(fit2 ,type = "response")
table(MyData$Survived,prob.prediksi>0.5)
##
## FALSE TRUE
## 0 389 35
## 1 112 176
dengan akurasi model 0,7935 atau 79%
Setelah itu kita melakukan uji wald untuk melihat parameter mana saja yang memberikan pengaruh yang signifikan terhadap variabel dependent.
summary(fit2)
##
## Call:
## glm(formula = Survived ~ Age + Pclass + Sex + Age * Pclass +
## Sex * Pclass, family = binomial, data = MyData)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2671 -0.6893 -0.4587 0.4873 2.3696
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.400099 0.725725 7.441 9.99e-14 ***
## Age -0.034001 0.012703 -2.677 0.00744 **
## Pclass -2.293156 0.408456 -5.614 1.97e-08 ***
## Sexmale -4.645139 0.604169 -7.688 1.49e-14 ***
## Age:Pclass -0.004016 0.008990 -0.447 0.65508
## Pclass:Sexmale 1.488765 0.353631 4.210 2.55e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 960.90 on 711 degrees of freedom
## Residual deviance: 623.67 on 706 degrees of freedom
## AIC: 635.67
##
## Number of Fisher Scoring iterations: 6
Berdasarkan output diatas, didapatkan hasil yang signifikan yaitu faktor Age, Pclass, Sexmale, Interaksi Pclass dan Sexmale.
Model regresi logistik parameter signifikan akan diinterpretasi menggunakan nilai oods ratio. Nilai odds ratio diberikan pada tabel berikut.
| Variabel Prediktor | Koefisien (β) | OR (Exp(β)) |
|---|---|---|
| Age | −0.034001 | 0,966570 |
| Pclass | −2,293156 | 0,100947 |
| Sexmale | −4,645139 | 9,6081 x 10−3 |
| Pclass*Sexmale | 1,488765 | 4,431619 |
Nilai odds ratio untuk umur sebesar 0,966570. Hal ini menunjukan bahwa odds bertahan hidup akan berkurang sebesar 0,966570 untuk setiap 1 satuan umur bertambah. Odds bertahan hidup untuk Passanger Class akan berkurang sebesar 0,100947 untuk setiap 1 kategori kelas bertambah. Odds bertahan hidup untuk jenis kelamin laki-laki adalah 9,6081 x 10−3 kali odds untuk jenis kelamin perempuan. Odds ratio untuk interaksi Passanger Class dan jenis kelamin laki-laki (Pclass*Sexmale) menunjukan untuk setiap kenaikan 1 kategori kelas pada jenis kelamin laki-laki akan memiliki odds bertahan hidup 4,431619 kali dibandingkan dengan jenis kelamin perempuan.
Berdasarkan hasil pengamatan yang sudah dilakukan berdasarkan materi analisis data kategorik maka dapat disimpulkan bahwa :