Pengenalan

Analisis ini berkaitan dengan mengeksplorasi status penyakit jantung dan memahami hubungannya dengan banyak variabel dan faktor yang berbeda. Secara khusus, bertujuan untuk pemahaman yang lebih baik tentang bagaimana jenis kelamin (pria dan wanita), usia, kadar kolesterol serum, nyeri dada, tekanan darah, dan detak jantung maksimum berperan dalam prevalensi penyakit jantung. Variabel tersebut akan dianalisis kemudian digunakan untuk membuat model Logistic Linear dan k-NN untuk memprediksi adanya kemungkinan menderita penyakit jantung.

Kita akan menggunakan dataset Heart Disease UCI dari Kaggle yang bisa dilihat disini

Apa yang akan kita lakukan ?

Kita akan mengolah data dengan Logistic Linear dan k-NN menggunakan dataset Heart Disease UCI. Langkah - langkah yang akan kita lakukan sebagai berikut:

  • Kita akan menetapkan target sebagai variabel target dan variabel lainnya sebagai prediktor
  • Kita cek korelasi variabel target dengan variabel prediktor lainnya
  • Membuat model Logistic Linear dan k-NN untuk menjelaskan variabel target berdasarkan prediktor
  • Menguji validitas model
  • Membuat evaluasi dari kedua model
  • Memberi kesimpulan dari kedua model

Persiapan Data

Import Library

library(class) #kNN
library(caret) #confusionMatrix
library(dplyr) #mengolah data sebelum di proses
library(DT) #datatables
library(GGally) #ggcorr
library(ggplot2)

Import Data

data <- read.csv("heart.csv")
glimpse(data)
## Rows: 303
## Columns: 14
## $ ï..age   <int> 63, 37, 41, 56, 57, 57, 56, 44, 52, 57, 54, 48, 49, 64, 58...
## $ sex      <int> 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0...
## $ cp       <int> 3, 2, 1, 1, 0, 0, 1, 1, 2, 2, 0, 2, 1, 3, 3, 2, 2, 3, 0, 3...
## $ trestbps <int> 145, 130, 130, 120, 120, 140, 140, 120, 172, 150, 140, 130...
## $ chol     <int> 233, 250, 204, 236, 354, 192, 294, 263, 199, 168, 239, 275...
## $ fbs      <int> 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0...
## $ restecg  <int> 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1...
## $ thalach  <int> 150, 187, 172, 178, 163, 148, 153, 173, 162, 174, 160, 139...
## $ exang    <int> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0...
## $ 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    <int> 0, 0, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 1, 2, 0, 2, 2...
## $ ca       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2...
## $ thal     <int> 1, 2, 2, 2, 2, 1, 2, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2...
## $ target   <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...

Dataset Heart Disease UCI memiliki variabel-variabel sebagai berikut:

  • ca = number of major vessels (0-3) colored by flourosopy
  • cp = chest pain type (0 - 3)
  • exang = exercise induced angina (1 = yes; 0 = no)
  • fbs = fbs(fasting blood sugar > 120 mg/dl) (1 = true; 0 = false)
  • restecg = resting electrocardiographic results
  • sex = sex(1 = male; 0 = female)
  • slope = the slope of the peak exercise ST segment
  • thal = 3 = normal; 6 = fixed defect; 7 = reversable defect
  • age = age in years
  • chol = cholserum cholestoral in mg/dl
  • oldpeak = oldpeakST depression induced by exercise relative to rest
  • testbps = resting blood pressure (in mm/Hg on admission to the hospital)
  • thalach = maximum heart rate achieved
  • target = target 1 or 0

Dataset memiliki 303 baris dan 14 kolom. Variabel target kita adalah target dan kita akan menggunakan variabel lain sebagai prediktor.

Data Wrangling

names(data)
##  [1] "ï..age"   "sex"      "cp"       "trestbps" "chol"     "fbs"     
##  [7] "restecg"  "thalach"  "exang"    "oldpeak"  "slope"    "ca"      
## [13] "thal"     "target"

Pada dataset terdapat nama kolom ï..age, kolom yang terdiri dari data umur target memiliki nama yang kurang rapih. Oleh karena itu mari kita rapikan terlebih dahulu dengan mengubah nama kolom menjadi age

colnames(data)[colnames(data)=="ï..age"] = "age"

Langkah pertama kita cek apakah terdapat missing value pada dateset

colSums(is.na(data))
##      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        0

Visualisasi Data

Kita akan memvisualisasikan variabel target dengan beberapa prediktornya

Kolerasi Antar Variabel

Lihat Kolerasi Antar Variabel

ggcorr(data = data, label = T, hjust = 0.9, cex=3)

Mari kita ubah tipe data variabel target menjadi factor

data$target <- ifelse(data$target == 1, "yes", "no")
data$target <- as.factor(data$target)

Variabel target dengan age

ggplot(data, mapping = aes( x = target, y = age, fill=target)) +
  geom_boxplot() +
  labs(x = "Target Status", 
       y = "Age",
       title = "Kolerasi Age dengan Target",
       subtitle = "Heart Disease UCI") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        legend.position = "none")

Variabel target dengan chol

ggplot(data, mapping = aes( x = target, y = chol, fill = target)) +
  geom_boxplot() +
  labs(x = "Target Status", 
       y = "Cholesterol Level (mg/dL)",
       title = "Cholesterol Levels Despite Heart Disease",
       subtitle = "Heart Disease UCI") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        legend.position = "none") +
   scale_fill_manual(values=c( "#E69F00", "#56B4E9"))

Analisis Data

  1. Cek Proporsi Target
prop.table(table(data$target))
## 
##        no       yes 
## 0.4554455 0.5445545

Proporsi kelas target seimbang.

  1. Cek range nilai antar prediktor (untuk k-NN)
summary(data)
##       age             sex               cp           trestbps    
##  Min.   :29.00   Min.   :0.0000   Min.   :0.000   Min.   : 94.0  
##  1st Qu.:47.50   1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:120.0  
##  Median :55.00   Median :1.0000   Median :1.000   Median :130.0  
##  Mean   :54.37   Mean   :0.6832   Mean   :0.967   Mean   :131.6  
##  3rd Qu.:61.00   3rd Qu.:1.0000   3rd Qu.:2.000   3rd Qu.:140.0  
##  Max.   :77.00   Max.   :1.0000   Max.   :3.000   Max.   :200.0  
##       chol            fbs            restecg          thalach     
##  Min.   :126.0   Min.   :0.0000   Min.   :0.0000   Min.   : 71.0  
##  1st Qu.:211.0   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:133.5  
##  Median :240.0   Median :0.0000   Median :1.0000   Median :153.0  
##  Mean   :246.3   Mean   :0.1485   Mean   :0.5281   Mean   :149.6  
##  3rd Qu.:274.5   3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:166.0  
##  Max.   :564.0   Max.   :1.0000   Max.   :2.0000   Max.   :202.0  
##      exang           oldpeak         slope             ca        
##  Min.   :0.0000   Min.   :0.00   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.00   1st Qu.:1.000   1st Qu.:0.0000  
##  Median :0.0000   Median :0.80   Median :1.000   Median :0.0000  
##  Mean   :0.3267   Mean   :1.04   Mean   :1.399   Mean   :0.7294  
##  3rd Qu.:1.0000   3rd Qu.:1.60   3rd Qu.:2.000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :6.20   Max.   :2.000   Max.   :4.0000  
##       thal       target   
##  Min.   :0.000   no :138  
##  1st Qu.:2.000   yes:165  
##  Median :2.000            
##  Mean   :2.314            
##  3rd Qu.:3.000            
##  Max.   :3.000

Modelling

Pisahkan data

Sebelumnya kita perlu membagi data menjadi dataset train dan dataset uji. Kita akan menggunakan dataset train untuk melatih model logistic linier. Dataset pengujian akan digunakan sebagai pembanding dan melihat apakah model menjadi overfit atau tidak dalam memprediksi data baru yang belum terlihat selama fase pelatihan.

Kita akan membagi sebanyak 80% dari data sebagai data pelatihan dan 20% sebagai data pengujian.

RNGkind(sample.kind = "Rounding")
set.seed(100)

intrain <- sample(nrow(data), nrow(data)*0.8)
data_train <- data[intrain,]
data_test <- data[-intrain,]

Kita cek kembali proporsi target pada data train

prop.table(table(data_train$target))
## 
##        no       yes 
## 0.4421488 0.5578512

proporsi kelas masih seimbang

Logistic Regression

model_log <- glm(target ~ ., data = data_train, family = "binomial")
summary(model_log)
## 
## Call:
## glm(formula = target ~ ., family = "binomial", data = data_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8007  -0.3274   0.1239   0.5334   2.4415  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.484818   3.088066   1.776 0.075711 .  
## age         -0.012713   0.026900  -0.473 0.636494    
## sex         -1.943239   0.559747  -3.472 0.000517 ***
## cp           1.177436   0.236588   4.977 6.47e-07 ***
## trestbps    -0.021032   0.012370  -1.700 0.089082 .  
## chol        -0.001761   0.004492  -0.392 0.695113    
## fbs          0.210562   0.651905   0.323 0.746699    
## restecg      0.474580   0.410174   1.157 0.247265    
## thalach      0.009584   0.012107   0.792 0.428579    
## exang       -1.220145   0.496660  -2.457 0.014022 *  
## oldpeak     -0.893392   0.306090  -2.919 0.003515 ** 
## slope        0.492960   0.472003   1.044 0.296300    
## ca          -0.653665   0.231756  -2.820 0.004795 ** 
## thal        -0.807547   0.338266  -2.387 0.016972 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 332.24  on 241  degrees of freedom
## Residual deviance: 154.28  on 228  degrees of freedom
## AIC: 182.28
## 
## Number of Fisher Scoring iterations: 6

Pada model tidak terindikasi adanya perfect separation

Mari kita coba interpretasikan salah satu prediktor pada model. Kali ini kita akan menginterpretasikan prediktor age

exp(-0.012713)
## [1] 0.9873675

Ini bearti setiap kenaikan 1 tahun umur menyebabkan 0.98 kali kemungkinan menderita penyakit jantung.

k-NN

Scaling menggunakan: z-score standardization

Pertama-tama, untuk k-NN harus dipisahkan antara data prediktor (x) & data label (y). Kemudian data prediktor akan discaling untuk data train maupun test (berdasarkan mean dan standar deviasi dari data train).

data_train_x <- data_train %>% select(-target)

data_test_x <- data_test %>% select(-target)

data_train_y <- data_train %>% select(target)
  
data_test_y <- data_test %>% select(target)
data_train_xs <- scale(x = data_train_x)
data_test_xs <- scale(x = data_test_x, 
                     center = attr(data_train_xs, "scaled:center"), 
                     scale = attr(data_train_xs, "scaled:scale"))

Selanjutnya kita cari nilai optimum k

round(sqrt(nrow(data_train)))
## [1] 16

Karena jumlah kelas target ada 2, maka kita ambil nilai 17 untuk dijadikan nilai optimum k

knn.Label <- knn(train = data_train_xs, 
                 test = data_test_xs,
                 cl = data_train_y$target,
                 k = 17)

head(knn.Label)
## [1] yes yes yes yes yes yes
## Levels: no yes

Prediksi

Logistic Regression

# cari probability
log.Risk <- predict(model_log, newdata = data_test, type = "response")

# tentukan kelas
log.Label <- ifelse(log.Risk > 0.5, "yes", "no")

# ubah label ke tipe faktor
log.Label <- as.factor(log.Label)

head(log.Label)
##   1   2   5  10  11  12 
## yes  no yes yes  no yes 
## Levels: no yes

Evaluasi Model

Logistic Regression

cm_log <- confusionMatrix(data = log.Label,
                reference = data_test_y$target,
                positive = "yes")
cm_log
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction no yes
##        no  23   7
##        yes  8  23
##                                           
##                Accuracy : 0.7541          
##                  95% CI : (0.6271, 0.8554)
##     No Information Rate : 0.5082          
##     P-Value [Acc > NIR] : 7.39e-05        
##                                           
##                   Kappa : 0.5083          
##                                           
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 0.7667          
##             Specificity : 0.7419          
##          Pos Pred Value : 0.7419          
##          Neg Pred Value : 0.7667          
##              Prevalence : 0.4918          
##          Detection Rate : 0.3770          
##    Detection Prevalence : 0.5082          
##       Balanced Accuracy : 0.7543          
##                                           
##        'Positive' Class : yes             
## 

kNN

cm_knn <- confusionMatrix(data = knn.Label, 
                reference = data_test$target,
                positive = "yes") 
cm_knn
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction no yes
##        no  22   3
##        yes  9  27
##                                          
##                Accuracy : 0.8033         
##                  95% CI : (0.6816, 0.894)
##     No Information Rate : 0.5082         
##     P-Value [Acc > NIR] : 1.809e-06      
##                                          
##                   Kappa : 0.6077         
##                                          
##  Mcnemar's Test P-Value : 0.1489         
##                                          
##             Sensitivity : 0.9000         
##             Specificity : 0.7097         
##          Pos Pred Value : 0.7500         
##          Neg Pred Value : 0.8800         
##              Prevalence : 0.4918         
##          Detection Rate : 0.4426         
##    Detection Prevalence : 0.5902         
##       Balanced Accuracy : 0.8048         
##                                          
##        'Positive' Class : yes            
## 

Kesimpulan

Pada kasus ini kita akan memfokuskan pada metrics recall / Sensitivity

Berikut ini hasil performa pada uji coba data test di metrics recall/Sensitivity untuk setiap model :

  • log : 0.76
  • kNN : 0.90

Berdasarkan hasil performa di data test, model k-NN yang menghasilkan metrics evaluasi dengan nilai paling baik. Namun dari model tersebut kita tidak bisa langsung menginterpretasinya. Jika ingin hasil yang mudah diinterpretasi maka kita akan memilih model logistic regression