Prediksi Pasien Penyakit Jantung
Objective
Pada kesempatan kali ini, saya akan mencoba melakukan prediksi terhadap pasien penyakit jantung pada suatu rumah sakit akan diprediksi sakit atau tidak berdasarkan kategori dari beberapa variabel penunjangnya. Algoritma yang akan saya gunakan yaitu menggunakan logistik regression dan k-nearest neighbor yang termasuk dalam supervised learning.
Library and Setup
Sebelum itu kita harus melakukan install.package() pada package dplyr, MASS, gtools, gmodels, class, ggplot2 pada R Studio. Apabila telah ter-install, maka lakukan pengaktifan package menggunakan library().
Logistic Regression
Data Import
Dataset yang akan saya gunakan yaitu data mengenai pasien yang terkena penyakit jantung berdasarkan beberapa karakteristik yang menyertai yang dapat Anda unduh lansung pada Kaggle.
## 'data.frame': 303 obs. of 14 variables:
## $ ï..age : int 63 37 41 56 57 57 56 44 52 57 ...
## $ sex : int 1 1 0 1 0 1 0 1 1 1 ...
## $ cp : int 3 2 1 1 0 0 1 1 2 2 ...
## $ trestbps: int 145 130 130 120 120 140 140 120 172 150 ...
## $ chol : int 233 250 204 236 354 192 294 263 199 168 ...
## $ fbs : int 1 0 0 0 0 0 0 0 1 0 ...
## $ restecg : int 0 1 0 1 1 1 0 1 1 1 ...
## $ thalach : int 150 187 172 178 163 148 153 173 162 174 ...
## $ exang : int 0 0 0 0 1 0 0 0 0 0 ...
## $ oldpeak : num 2.3 3.5 1.4 0.8 0.6 0.4 1.3 0 0.5 1.6 ...
## $ slope : int 0 0 2 2 2 1 1 2 2 2 ...
## $ ca : int 0 0 0 0 0 0 0 0 0 0 ...
## $ thal : int 1 2 2 2 2 1 2 3 3 2 ...
## $ target : int 1 1 1 1 1 1 1 1 1 1 ...Informasi penting dalam data :
ï..age : dalam beberapa tahun
sex : (1 = laki-laki; 0 = perempuan)
cp : tipe nyeri yang paling parah
trestbps : melacak tekanan darah(dalam mm Hg saat masuk ke rumah sakit)
chol : kolestoral dalam mg / dl
fbs : (gula darah puasa> 120 mg / dl) (1 = benar; 0 = salah)
restecg : mengembalikan hasil elektrokardiografi
thalach : denyut jantung maksimum tercapai
exang : exercise induced angina (1 = ya; 0 = tidak)
oldpeak : ST depresi yang disebabkan oleh olahraga relatif terhadap istirahat
slope : kemiringan segmen ST latihan puncak
ca : jumlah pembuluh darah utama (0-3) diwarnai dengan fluoroskopi
thal : 3 = normal; 6 = cacat tetap; 7 = cacat yang dapat dibalik
target : 1 = sakit atau 0 = tidak sakit
Berikut ini gambaran sedikit pada data yang digunakan.
Data Manipulation
Pada beberapa variabel yang digunakan, terdapat ketidak sesuaian tipe data, oleh karena itu yang perlu kita lakukan adalah melakukan penyesuaian tipe data pada beberapa variabel yang ada.
heart <- heart %>%
mutate_if(is.integer, as.factor) %>%
mutate(sex = factor(sex, levels = c(0,1), labels = c("Female", "Male")),
fbs =factor(fbs, levels = c(0,1), labels = c("False", "True")),
exang = factor(exang, levels = c(0,1), labels = c("No", "Yes")),
target = factor(target, levels = c(0,1),
labels = c("Health", "Not Health")))
glimpse(heart)## Observations: 303
## Variables: 14
## $ ï..age <fct> 63, 37, 41, 56, 57, 57, 56, 44, 52, 57, 54, 48, 49, 64, 58...
## $ sex <fct> Male, Male, Female, Male, Female, Male, Female, Male, Male...
## $ cp <fct> 3, 2, 1, 1, 0, 0, 1, 1, 2, 2, 0, 2, 1, 3, 3, 2, 2, 3, 0, 3...
## $ trestbps <fct> 145, 130, 130, 120, 120, 140, 140, 120, 172, 150, 140, 130...
## $ chol <fct> 233, 250, 204, 236, 354, 192, 294, 263, 199, 168, 239, 275...
## $ fbs <fct> True, False, False, False, False, False, False, False, Tru...
## $ restecg <fct> 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1...
## $ thalach <fct> 150, 187, 172, 178, 163, 148, 153, 173, 162, 174, 160, 139...
## $ exang <fct> No, No, No, No, Yes, No, No, No, No, No, No, No, No, Yes, ...
## $ oldpeak <dbl> 2.3, 3.5, 1.4, 0.8, 0.6, 0.4, 1.3, 0.0, 0.5, 1.6, 1.2, 0.2...
## $ slope <fct> 0, 0, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 1, 2, 0, 2, 2...
## $ ca <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2...
## $ thal <fct> 1, 2, 2, 2, 2, 1, 2, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2...
## $ target <fct> Not Health, Not Health, Not Health, Not Health, Not Health...Selanjutnya yaitu melakukan pengecekan terhadap missing value. Missing value perlu kita cek terlebih dahulu agar tidak mengganggu dalam melakukan pemodelan nantinya.
## ï..age sex cp trestbps chol fbs restecg thalach
## 0 0 0 0 0 0 0 0
## exang oldpeak slope ca thal target
## 0 0 0 0 0 0Pre-Processing Data
Sebelum melakukan pemodelan, kita perlu melihat terlebih dahulu proporsi dari target variabel yang kita miliki pada kolom target.
##
## Health Not Health
## 0.4554455 0.5445545##
## Health Not Health
## 138 165Jika dilihat dari proporsi kedua kelas, sudah cukup seimbang, sehingga kita tidak terlalu membutuhkan pre-processing tambahan untuk menyeimbangkan proporsi antar dua kelas target variabel.
Splitting Train-Test
Langkah selanjutnya yaitu melakukan splitting train test data. Tujuannya yaitu pada data train akan kita gunakan untuk modeling, sedangkan data test akan kita gunakan sebagai penguji model yang sudah kita buat jika dihadapkan dengan unseen data. Selain itu hal ini dapat digunakan untuk melihat kemampuan model yang kita buat dalam menghadapi unseen data.
set.seed(303)
intrain <- sample(nrow(heart), nrow(heart)*0.7)
heart_train <- heart[intrain,]
heart_test <- heart[-intrain,]
heart$target %>%
levels()## [1] "Health" "Not Health"Modelling
Melakukan pemodelan menggunakan regresi logistik. Pemodelan menggunakan fungsi glm() dalam memodelkan menggunakan regresi logistik. Variabel yang digunakan adalah beberapa variabel yang kita anggap mempengaruhi target variabel, dimana variabel target menjadi variabel responnya.
model <- glm(formula = target~sex+cp+fbs+exang+oldpeak+slope+ca+thal, family = "binomial",
data = heart_train)
summary(model)##
## Call:
## glm(formula = target ~ sex + cp + fbs + exang + oldpeak + slope +
## ca + thal, family = "binomial", data = heart_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9452 -0.3315 0.1082 0.5169 2.9240
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2003 3.4150 -0.059 0.95322
## sexMale -0.9320 0.5862 -1.590 0.11184
## cp1 0.4605 0.6375 0.722 0.47004
## cp2 1.9331 0.6323 3.057 0.00223 **
## cp3 1.8515 0.7319 2.530 0.01142 *
## fbsTrue 0.4894 0.6819 0.718 0.47296
## exangYes -1.2490 0.5154 -2.423 0.01538 *
## oldpeak -0.3834 0.2632 -1.457 0.14514
## slope1 -0.6466 0.9489 -0.681 0.49559
## slope2 1.1892 1.0309 1.154 0.24865
## ca1 -1.9168 0.5763 -3.326 0.00088 ***
## ca2 -2.8735 0.8966 -3.205 0.00135 **
## ca3 -1.9388 1.1738 -1.652 0.09859 .
## ca4 15.6796 1072.3805 0.015 0.98833
## thal1 1.4827 3.3717 0.440 0.66011
## thal2 2.3340 3.2627 0.715 0.47439
## thal3 0.7604 3.2621 0.233 0.81568
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 292.00 on 211 degrees of freedom
## Residual deviance: 134.57 on 195 degrees of freedom
## AIC: 168.57
##
## Number of Fisher Scoring iterations: 15Model Fitting
Pada pemodelan yang pertama, masih banyak variabel prediktor yang tidak signifikan terhadap target variabel, oleh karena itu kita akan coba melakukan model fitting menggunakan metode stepwise.
## Start: AIC=168.57
## target ~ sex + cp + fbs + exang + oldpeak + slope + ca + thal
##
## Df Deviance AIC
## - fbs 1 135.09 167.09
## <none> 134.57 168.57
## - oldpeak 1 136.80 168.80
## - sex 1 137.13 169.13
## - exang 1 140.50 172.50
## - thal 3 144.57 172.57
## - cp 3 148.42 176.42
## - slope 2 146.97 176.97
## - ca 4 157.84 183.84
##
## Step: AIC=167.09
## target ~ sex + cp + exang + oldpeak + slope + ca + thal
##
## Df Deviance AIC
## <none> 135.09 167.09
## - oldpeak 1 137.41 167.41
## - sex 1 137.62 167.62
## - exang 1 140.77 170.77
## - thal 3 144.84 170.84
## - slope 2 147.52 175.52
## - cp 3 150.66 176.66
## - ca 4 157.84 181.84Dengan menggunakan metode backward pada stepwise, kita memperoleh model sebagai berikut.
##
## Call:
## glm(formula = target ~ sex + cp + exang + oldpeak + slope + ca +
## thal, family = "binomial", data = heart_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9721 -0.3359 0.1129 0.5304 2.8961
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.0353 3.8113 -0.009 0.992611
## sexMale -0.9220 0.5840 -1.579 0.114386
## cp1 0.5168 0.6277 0.823 0.410366
## cp2 2.0309 0.6173 3.290 0.001003 **
## cp3 1.9043 0.7299 2.609 0.009080 **
## exangYes -1.2164 0.5123 -2.375 0.017569 *
## oldpeak -0.3903 0.2632 -1.483 0.138067
## slope1 -0.6305 0.9560 -0.660 0.509564
## slope2 1.1990 1.0407 1.152 0.249272
## ca1 -1.8824 0.5687 -3.310 0.000934 ***
## ca2 -2.7818 0.8828 -3.151 0.001627 **
## ca3 -1.7568 1.1208 -1.567 0.117016
## ca4 15.6050 1070.2774 0.015 0.988367
## thal1 1.3575 3.7769 0.359 0.719276
## thal2 2.1320 3.6691 0.581 0.561197
## thal3 0.5864 3.6727 0.160 0.873156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 292.00 on 211 degrees of freedom
## Residual deviance: 135.09 on 196 degrees of freedom
## AIC: 167.09
##
## Number of Fisher Scoring iterations: 15Prediksi
Dengan menggunakan model2 hasil dari stepwise, kita akan coba prediksi menggunakan data test yang sudah kita miliki.
Kita akan coba melihat sebaran peluang prediksi data.
ggplot(heart_test, aes(x=prob_heart)) +
geom_density(lwd=0.5) +
labs(title = "Distribution of Probability Prediction Data") +
theme_minimal()
Pada grafik diatas, dapat diinterpretasikan bahwa hasil prediksi yang dilakukan lebih condong ke arah 1 yang artinya Not Health.
heart_test$pred_heart <- factor(ifelse(heart_test$prob_heart > 0.5, "Not Health","Health"))
heart_test[1:10, c("pred_heart", "target")]Dalam syntax diatas, ketika probabilitas data test lebih dari 0.5, artinya dia sakit atau Not Health.
Model Evaluation
Untuk mengevaluasi model yang telah kita buat, kita akan menggunakan confusion matrix.
library(caret)
log_conf <- confusionMatrix(heart_test$pred_heart, heart_test$target, positive = "Not Health")
log_conf## Confusion Matrix and Statistics
##
## Reference
## Prediction Health Not Health
## Health 33 4
## Not Health 9 45
##
## Accuracy : 0.8571
## 95% CI : (0.7681, 0.9217)
## No Information Rate : 0.5385
## P-Value [Acc > NIR] : 1.091e-10
##
## Kappa : 0.7101
##
## Mcnemar's Test P-Value : 0.2673
##
## Sensitivity : 0.9184
## Specificity : 0.7857
## Pos Pred Value : 0.8333
## Neg Pred Value : 0.8919
## Prevalence : 0.5385
## Detection Rate : 0.4945
## Detection Prevalence : 0.5934
## Balanced Accuracy : 0.8520
##
## 'Positive' Class : Not Health
## - Re-call/Sensitivity = dari semua data aktual yang positif, seberapa mampu proporsi model saya menebak benar.
- Specificity = dari semua data aktual yang negatif, seberapa mampu proporsi model saya menebak yang benar.
- Accuracy = seberapa mampu model saya menebak dengan benar target Y.
- Precision = dari semua hasil prediksi, seberapa mampu model saya dapat menebak benar kelas positif.
Re-call = (49)/(7+49) Specificity = (27)/(27+8) Accuracy = (27+49)/(27+7+8+49) Precision = (49)/(49+8)
Recall <- round((27)/(27+8),2)
Specificity <- round((49)/(7+49),2)
Accuracy <- round((27+49)/(27+7+8+49),2)
Precision <- round((49)/(49+8),2)
performance <- cbind.data.frame(Accuracy, Recall, Precision, Specificity)
performanceBerdasarkan hasil confusionMatrix diatas, dapat kita ambil informasi bahwa kemampuan model dalam menebak target Y (Health dan Not Health) sebesar 83,5%. Sedangkan dari keluruhan data aktual orang yang Not Health, model dapat mampu menebak benar sebesar 87,5%. Dari keseluruhan data aktual orang yang Health, model mampu menebak dengan benar sebesar 77,1%. Dari keseluruhan hasil prediksi yang mampu ditebak oleh model, model mampu menebak benar kelas positif sebesar 86%.
Tuning Cutoff
Digunakan untuk mengetahui threshold maksimum dari apa yang akan kita teliti.
# tuning cutoff
performa <- function(cutoff, prob, ref, postarget, negtarget)
{
predict <- factor(ifelse(prob >= cutoff, postarget, negtarget))
conf <- caret::confusionMatrix(predict , ref, positive = postarget)
acc <- conf$overall[1]
rec <- conf$byClass[1]
prec <- conf$byClass[3]
spec <- conf$byClass[2]
mat <- t(as.matrix(c(rec , acc , prec, spec)))
colnames(mat) <- c("recall", "accuracy", "precicion", "specificity")
return(mat)
}
co <- seq(0.01,0.80,length=100)
result <- matrix(0,100,4)
for(i in 1:100){
result[i,] = performa(cutoff = co[i],
prob = heart_test$prob_heart,
ref = heart_test$target,
postarget = "Not Health",
negtarget = "Health")
}
data_frame("Recall" = result[,1],
"Accuracy" = result[,2],
"Precision" = result[,3],
"Specificity" = result[,4],
"Cutoff" = co) %>%
gather(key = "performa", value = "value", 1:4) %>%
ggplot(aes(x = Cutoff, y = value, col = performa)) +
geom_line(lwd = 1.5) +
scale_color_manual(values = c("darkred","darkgreen","orange", "blue")) +
scale_y_continuous(breaks = seq(0,1,0.1), limits = c(0,1)) +
scale_x_continuous(breaks = seq(0,1,0.1)) +
labs(title = "Tradeoff model perfomance") +
theme_minimal() +
theme(legend.position = "top",
panel.grid.minor.y = element_blank(),
panel.grid.minor.x = element_blank())
Berdasarkan Tradeoff model performance diatas, dapat kita tahu bahwa dengan cutoff 0.5 kita memperoleh nilai precision dan re-call yang agak tinggi, namun nilai accuracy dan nilai specificity agak rendah.
Model Interpretation
Interpretasi model : 1. Odds orang berjenis kelamin laki-laki = 0,206 < 1. Artinya orang dengan jenis kelamin laki-laki lebih kecil 20,6% untuk Not Health dibandingkan dengan perempuan.
K-Nearest Neighbour
Pre-Processing Data
Membuat variabel dummy dari data-data kategori yang akan digunakan dalam klasifikasi.
dmy <- dummyVars(" ~target+sex+cp+fbs+exang+oldpeak+slope+ca+thal", data = heart)
dmy <- data.frame(predict(dmy, newdata = heart))
str(dmy)## 'data.frame': 303 obs. of 25 variables:
## $ target.Health : num 0 0 0 0 0 0 0 0 0 0 ...
## $ target.Not.Health: num 1 1 1 1 1 1 1 1 1 1 ...
## $ sex.Female : num 0 0 1 0 1 0 1 0 0 0 ...
## $ sex.Male : num 1 1 0 1 0 1 0 1 1 1 ...
## $ cp.0 : num 0 0 0 0 1 1 0 0 0 0 ...
## $ cp.1 : num 0 0 1 1 0 0 1 1 0 0 ...
## $ cp.2 : num 0 1 0 0 0 0 0 0 1 1 ...
## $ cp.3 : num 1 0 0 0 0 0 0 0 0 0 ...
## $ fbs.False : num 0 1 1 1 1 1 1 1 0 1 ...
## $ fbs.True : num 1 0 0 0 0 0 0 0 1 0 ...
## $ exang.No : num 1 1 1 1 0 1 1 1 1 1 ...
## $ exang.Yes : num 0 0 0 0 1 0 0 0 0 0 ...
## $ oldpeak : num 2.3 3.5 1.4 0.8 0.6 0.4 1.3 0 0.5 1.6 ...
## $ slope.0 : num 1 1 0 0 0 0 0 0 0 0 ...
## $ slope.1 : num 0 0 0 0 0 1 1 0 0 0 ...
## $ slope.2 : num 0 0 1 1 1 0 0 1 1 1 ...
## $ ca.0 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ ca.1 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ca.2 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ca.3 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ca.4 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ thal.0 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ thal.1 : num 1 0 0 0 0 1 0 0 0 0 ...
## $ thal.2 : num 0 1 1 1 1 0 1 0 0 1 ...
## $ thal.3 : num 0 0 0 0 0 0 0 1 1 0 ...Menghapus variabel dummy yang variabel sebelumnya hanya terdapat 2 kategori.
Mengetahui nama-nama dari variabel dummy yang terbentuk.
## [1] "target.Not.Health" "sex.Male" "cp.0"
## [4] "cp.1" "cp.2" "cp.3"
## [7] "fbs.True" "exang.Yes" "oldpeak"
## [10] "slope.0" "slope.1" "slope.2"
## [13] "ca.0" "ca.1" "ca.2"
## [16] "ca.3" "ca.4" "thal.0"
## [19] "thal.1" "thal.2" "thal.3"Membentuk data training dan data testing dari data dmy yang telah terbentuk.
set.seed(300)
dmy_train <- dmy[intrain,2:21]
dmy_test <- dmy[-intrain,2:21]
dmy_train_label <- dmy[intrain,1]
dmy_test_label <- dmy[-intrain,1]Melakukan prediksi dengan K-NN
Membuat confusion matriks dari prediski K-NN
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 33 5
## 1 9 44
##
## Accuracy : 0.8462
## 95% CI : (0.7554, 0.9133)
## No Information Rate : 0.5385
## P-Value [Acc > NIR] : 5.323e-10
##
## Kappa : 0.6884
##
## Mcnemar's Test P-Value : 0.4227
##
## Sensitivity : 0.8980
## Specificity : 0.7857
## Pos Pred Value : 0.8302
## Neg Pred Value : 0.8684
## Prevalence : 0.5385
## Detection Rate : 0.4835
## Detection Prevalence : 0.5824
## Balanced Accuracy : 0.8418
##
## 'Positive' Class : 1
## Berdasarkan hasil confusion matrix diatas, dapat kita ketahui bahwa kemampuan model dalam menebak target Y sebesar 87,9%. Sedangkan berdasarkan data aktual orang yang memiliki status Not Health, model dapat menebak dengan benar sebesar 94,6%. berdasarkan data aktual orang yang memiliki status Health, model dapat menebak dengan benar sebesar 77,14%. Dari keseluruhan hasil prediksi yang mampu ditebak oleh model, model mampu menebak benar kelas positif sebesar 86,9%.
Model Evaluation Logistic Regression and K-NN
eval_logit <- data_frame(Accuracy = log_conf$overall[1],
Recall = log_conf$byClass[1],
Specificity = log_conf$byClass[2],
Precision = log_conf$byClass[3])
eval_knn <- data_frame(Accuracy = pred_knn_conf$overall[1],
Recall = pred_knn_conf$byClass[1],
Specificity = log_conf$byClass[2],
Precision = pred_knn_conf$byClass[3])Jika dilihat dari kedua metode tersebut, yaitu dengan menggunakan Regresi Logistik dan K-NN, kemampupuan model dalam memprediksi benar dari data aktual orang yang Not Health lebih baik dengan menggunakan metode K-NN karena memiliki nilai precision= 86,9% lebih besar dari pada menggunakan metode regresi logistik.
Conclusion
Jika saya mengibaratkan diri saya seorang dokter penyakit jantung, dimana treatment yang akan saya lakukan ke pasien saya yang sakit jantung dengan yang tidak sakit jantung sangat berbeda. Treatment yang akan diberikan untuk pasien yang sakit jantung jika diberikan kepada pasien yang tidak sakit, akan membahayakan pasien yang tidak sakit tersebut. Oleh karena itu, saya akan sangat melihat metric precision yang ada, dimana saya tidak ingin model yang saya buat salah dalam memprediksi mana pasien yang benar-benar sakit jantung atau yang tidak sakit jantung.