1 Probabilitas Pengidap Penyakit Jantung Menggunakan Logistic Regression

disusun oleh Reza Lutfi Ismail

knitr::include_graphics("hd.jfif")

1.1 Background

Project kali ini menjelaskan tentang prediksi pasien yang memiliki kemungkinan untuk mengidap penyakit jantung berdasarkan data. adapun, data yang didapatkan adalah data dari kaggle pada link: https://www.kaggle.com/ronitf/heart-disease-uci

2 Data Preprocessing

2.1 Import Library

library(tidyverse) # data processing 
library(gtools)
library(class) # classification

2.2 Import Data

heart <- read_csv("data_input/heart.csv")
head(heart)
#> # A tibble: 6 x 14
#>     age   sex    cp trestbps  chol   fbs restecg thalach exang oldpeak slope
#>   <dbl> <dbl> <dbl>    <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>   <dbl> <dbl>
#> 1    63     1     3      145   233     1       0     150     0     2.3     0
#> 2    37     1     2      130   250     0       1     187     0     3.5     0
#> 3    41     0     1      130   204     0       0     172     0     1.4     2
#> 4    56     1     1      120   236     0       1     178     0     0.8     2
#> 5    57     0     0      120   354     0       1     163     1     0.6     2
#> 6    57     1     0      140   192     0       1     148     0     0.4     1
#> # ... with 3 more variables: ca <dbl>, thal <dbl>, target <dbl>

2.3 EDA

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

Atribut Informasi:

  • age = usia
  • sex = gender (1 = laki-laki, 0 = perempuan)
  • cp = tipe nyeri paling parah (4 values)
  • trestbps = melacak tekanan resting blood pressure
  • chol = serum kolesterol dalam mg/dl
  • fbs = gula darah puasa > 120 mg/dl (1 = ya, 0 = tidak)
  • restecg = mengembalikan hasil electrocardiographic results (values 0,1,2)
  • thalach = maksimum denyut jantung yang tercapai
  • exang = latihan angina yang diinduksi (1 = ya, 0 = tidak)
  • oldpeak = ST depresi yang disebabkan olahraga relatif terhadap istirahat
  • slope = kemiringan segmen ST latihan puncak
  • ca = jumlah pembuluh darah besar (0-3) warna berdasarkan flourosopy
  • thal = (3 = normal; 6 = cacat selamanya; 7 = cacat yang dapat sembuh)
  • target = (1 = sakit; 0 = tidak sakit)

berdasarkan informasi dari type data, terdapat beberapa data yang berulang yaitu: - cp - slope - ca - thal - sex - fbs - exang - target

variabel tersebut nantinya akan dijadikan factor (kategorikal)

heart_clean <- heart %>% 
  select(-age) %>% # menghapus kolom age 
  mutate(cp = as.factor(cp), # ubah kolom cp menjdi kategorikal
         slope = as.factor(slope), # ubah kolom slope menjadi kategorikal
         ca = as.factor(ca), # ubah kolom ca menjadi kategorikal 
         thal = as.factor(thal)) %>% # ubah kolom thal menjadi kategorikal
  mutate(sex = factor(sex, levels = c(0,1), labels = c("perempuan", "laki-laki")), # ubah kolom sex menjadi kategorikal
         fbs = factor(fbs, levels = c(0,1), labels = c("false", "true")), # ubah kolom fbs menjadi kategorikal
         exang = factor(exang, levels = c(0,1), labels = c("no", "yes")), # ubah kolom exang menjadi kategorikal
         target = factor(target, levels = c(0,1), 
                         labels = c("health", "not health"))) # ubah kolom target menjadi kategorikal

Cek 6 baris data frame heart_clean

# check data frame heart_clean
head(heart_clean)
#> # A tibble: 6 x 13
#>   sex   cp    trestbps  chol fbs   restecg thalach exang oldpeak slope ca   
#>   <fct> <fct>    <dbl> <dbl> <fct>   <dbl>   <dbl> <fct>   <dbl> <fct> <fct>
#> 1 laki~ 3          145   233 true        0     150 no        2.3 0     0    
#> 2 laki~ 2          130   250 false       1     187 no        3.5 0     0    
#> 3 pere~ 1          130   204 false       0     172 no        1.4 2     0    
#> 4 laki~ 1          120   236 false       1     178 no        0.8 2     0    
#> 5 pere~ 0          120   354 false       1     163 yes       0.6 2     0    
#> 6 laki~ 0          140   192 false       1     148 no        0.4 1     0    
#> # ... with 2 more variables: thal <fct>, target <fct>

Cek perbandingan antara peluang target yang health dan not heatlh

prop.table(table(heart_clean$target))
#> 
#>     health not health 
#>  0.4554455  0.5445545

dari informasi diatas, diketahui bahwa perbandingan antara target health dan not health berimbang. artinya sebaran data sudah baik

2.4 Splitting Data Train dan Data Test

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

index <- sample(nrow(heart_clean), 0.8*nrow(heart_clean))
data_train <- heart_clean[index, ]
data_test <- heart_clean[-index, ]

3 Modelling

3.1 Model All Predictor

membuat model dengan menggunakan semua prediktor pada dataframe heart_clean

model_all <- glm(formula = target ~ sex + cp + fbs + exang + oldpeak + slope + ca + thal, family = "binomial", data = data_train)
summary(model_all)
#> 
#> Call:
#> glm(formula = target ~ sex + cp + fbs + exang + oldpeak + slope + 
#>     ca + thal, family = "binomial", data = data_train)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -2.7263  -0.3748   0.1274   0.4629   2.7648  
#> 
#> Coefficients:
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)   -12.9718  1455.3983  -0.009 0.992889    
#> sexlaki-laki   -1.4307     0.5361  -2.669 0.007615 ** 
#> cp1             1.2477     0.6199   2.013 0.044143 *  
#> cp2             2.2604     0.5893   3.836 0.000125 ***
#> cp3             2.0516     0.6959   2.948 0.003200 ** 
#> fbstrue         0.3536     0.6144   0.576 0.564908    
#> exangyes       -0.4509     0.4723  -0.955 0.339801    
#> oldpeak        -0.4456     0.2543  -1.752 0.079732 .  
#> slope1         -0.6207     0.9231  -0.672 0.501357    
#> slope2          0.7462     1.0176   0.733 0.463385    
#> ca1            -1.9731     0.5342  -3.693 0.000221 ***
#> ca2            -2.8602     0.8091  -3.535 0.000407 ***
#> ca3            -2.2190     0.8905  -2.492 0.012711 *  
#> ca4             0.8897     1.5436   0.576 0.564359    
#> thal1          14.9188  1455.3979   0.010 0.991821    
#> thal2          15.0878  1455.3978   0.010 0.991729    
#> thal3          13.7363  1455.3977   0.009 0.992470    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 333.83  on 241  degrees of freedom
#> Residual deviance: 159.89  on 225  degrees of freedom
#> AIC: 193.89
#> 
#> Number of Fisher Scoring iterations: 14

3.2 Backward

setelah dibuat model dengan menggunakan semua prediktor, kita akan melakukan step-wise backward untuk menentukan prediktor mana saja yang memiliki pengaruh besar terhadap model yang kita miliki, lalu di assign kedalam object backward

backward <- step(object = model_all, direction = "backward")
#> Start:  AIC=193.89
#> target ~ sex + cp + fbs + exang + oldpeak + slope + ca + thal
#> 
#>           Df Deviance    AIC
#> - fbs      1   160.22 192.22
#> - exang    1   160.79 192.79
#> <none>         159.89 193.89
#> - oldpeak  1   163.13 195.13
#> - slope    2   167.66 197.66
#> - thal     3   170.02 198.02
#> - sex      1   167.52 199.52
#> - cp       3   180.23 208.23
#> - ca       4   189.93 215.93
#> 
#> Step:  AIC=192.22
#> target ~ sex + cp + exang + oldpeak + slope + ca + thal
#> 
#>           Df Deviance    AIC
#> - exang    1   161.14 191.14
#> <none>         160.22 192.22
#> - oldpeak  1   163.80 193.80
#> - slope    2   167.80 195.80
#> - thal     3   170.30 196.30
#> - sex      1   167.94 197.94
#> - cp       3   181.94 207.94
#> - ca       4   190.05 214.05
#> 
#> Step:  AIC=191.14
#> target ~ sex + cp + oldpeak + slope + ca + thal
#> 
#>           Df Deviance    AIC
#> <none>         161.14 191.14
#> - oldpeak  1   165.28 193.28
#> - slope    2   169.26 195.26
#> - thal     3   172.38 196.38
#> - sex      1   168.62 196.62
#> - ca       4   190.79 212.79
#> - cp       3   191.42 215.42

didapat bahwa prediktor yang digunakan adalah: - sex - cp - oldpeak - slope - ca - thal

lalu, kita akan buat kolom baru untuk prediksi model backward kedalam dataframe data_test dan kita namakan sebagai pred_target

data_test$pred_target <- predict(object = backward, newdata = data_test, type = "response")

maka, didapat dataframe data_test menjadi sebagai berikut:

data_test[,c("target", "pred_target")]
#> # A tibble: 61 x 2
#>    target     pred_target
#>    <fct>            <dbl>
#>  1 not health       0.968
#>  2 not health       0.908
#>  3 not health       0.383
#>  4 not health       0.781
#>  5 not health       0.891
#>  6 not health       0.665
#>  7 not health       0.988
#>  8 not health       0.631
#>  9 not health       0.147
#> 10 not health       0.972
#> # ... with 51 more rows

untuk menentukan klasifikasi prediksi target yang health dan not health, kita tentukan treshold jika nilai pred_target > 0.5 maka target adalah not health

data_test$pred <- factor(ifelse(data_test$pred_target > 0.5, "not health", "health"))
data_test[1:10, c("pred","target")]
#> # A tibble: 10 x 2
#>    pred       target    
#>    <fct>      <fct>     
#>  1 not health not health
#>  2 not health not health
#>  3 health     not health
#>  4 not health not health
#>  5 not health not health
#>  6 not health not health
#>  7 not health not health
#>  8 not health not health
#>  9 health     not health
#> 10 not health not health

3.3 Evaluate Model

untuk mengetahui bahwa model yang kita memiliki performa yang baik atau tidak, maka akan kita lakukan evaluasi model dengan menggunakan confusionmatrix terhadap varibel target dan pred

library(caret)
conf <- confusionMatrix(data_test$pred, data_test$target, positive = "not health")
conf
#> Confusion Matrix and Statistics
#> 
#>             Reference
#> Prediction   health not health
#>   health         21          4
#>   not health      6         30
#>                                           
#>                Accuracy : 0.8361          
#>                  95% CI : (0.7191, 0.9185)
#>     No Information Rate : 0.5574          
#>     P-Value [Acc > NIR] : 0.000003844     
#>                                           
#>                   Kappa : 0.6652          
#>                                           
#>  Mcnemar's Test P-Value : 0.7518          
#>                                           
#>             Sensitivity : 0.8824          
#>             Specificity : 0.7778          
#>          Pos Pred Value : 0.8333          
#>          Neg Pred Value : 0.8400          
#>              Prevalence : 0.5574          
#>          Detection Rate : 0.4918          
#>    Detection Prevalence : 0.5902          
#>       Balanced Accuracy : 0.8301          
#>                                           
#>        'Positive' Class : not health      
#>

berdasarkan informasi diatas, bahwa kemampuan model untuk menebak target (health dan not health) sebesar 84%, sedangkan dari keseluruhan data aktual orang yang not health, model mampu menebak sebesar 78%. sedangkan dari keseluruhan data aktual orang yang health, model mampu menebak sebesar 88%. sedangkan dari data prediksi, model mampu menebak kelas positif (not health) sebesar 83%.

3.4 Model Interpretation

exp(model_all$coefficients) %>% 
  data.frame()
#>                                 .
#> (Intercept)        0.000002324976
#> sexlaki-laki       0.239134601803
#> cp1                3.482396105455
#> cp2                9.587034101791
#> cp3                7.780103846680
#> fbstrue            1.424193489568
#> exangyes           0.637080708374
#> oldpeak            0.640424036394
#> slope1             0.537589128457
#> slope2             2.108944004662
#> ca1                0.139032149025
#> ca2                0.057256991607
#> ca3                0.108722912959
#> ca4                2.434328494055
#> thal1        3013951.305251823738
#> thal2        3568935.718299917411
#> thal3         923820.034745599725

4 Conclusion

berdasarkan informasi dari model yang telah dibuat, didapat bahwa laki-laki memiliki probabilitas mengidap not health 24% lebih kecil dibandingkan dengan perempuan.