2502026873_FinalExam

Link Video Penjelasan : https://youtu.be/fKa1mHW8hNs

Load Data

library(readr)

## Warning: package 'readr' was built under R version 4.1.3

library(ggplot2) 

## Warning: package 'ggplot2' was built under R version 4.1.3

library(Hmisc)

## Warning: package 'Hmisc' was built under R version 4.1.3

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':
## 
##     format.pval, units

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.1.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ROCR)

## Warning: package 'ROCR' was built under R version 4.1.3

library(rpart)
library(rpart.plot)

## Warning: package 'rpart.plot' was built under R version 4.1.3

library(caret)

## Warning: package 'caret' was built under R version 4.1.3

## 
## Attaching package: 'caret'

## The following object is masked from 'package:survival':
## 
##     cluster

path <- "D:/R/heartData.csv"
heartdf <- read.csv(path)

1A. Exploratory Data Analysis

dim(heartdf)

## [1] 918  12

str(heartdf)

## 'data.frame':    918 obs. of  12 variables:
##  $ Age           : int  40 49 37 48 54 39 45 54 37 48 ...
##  $ Sex           : chr  "M" "F" "M" "F" ...
##  $ ChestPainType : chr  "ATA" "NAP" "ATA" "ASY" ...
##  $ RestingBP     : int  140 160 130 138 150 120 130 110 140 120 ...
##  $ Cholesterol   : int  289 180 283 214 195 339 237 208 207 284 ...
##  $ FastingBS     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ RestingECG    : chr  "Normal" "Normal" "ST" "Normal" ...
##  $ MaxHR         : int  172 156 98 108 122 170 170 142 130 120 ...
##  $ ExerciseAngina: chr  "N" "N" "N" "Y" ...
##  $ Oldpeak       : num  0 1 0 1.5 0 0 0 0 1.5 0 ...
##  $ ST_Slope      : chr  "Up" "Flat" "Up" "Flat" ...
##  $ HeartDisease  : int  0 1 0 1 0 0 0 0 1 0 ...

EXPLANATION:

Dataset yang saya gunakan memiliki 918 baris dan 12 atribut. Selanjutnya, saya melihat data teratas dari dataset yang saya gunakan. Dapat dilihat, ada berbagai atribut pada dataset terdiri dari tipe data integer, character, dan numeric.

heartdf$Sex <- as.factor(heartdf$Sex)
heartdf$ExerciseAngina <- as.factor(heartdf$ExerciseAngina)
heartdf$RestingECG <- as.factor(heartdf$RestingECG)
heartdf$ChestPainType <- as.factor(heartdf$ChestPainType)
heartdf$ST_Slope <- as.factor(heartdf$ST_Slope)
heartdf$FastingBS <- as.factor(heartdf$FastingBS)
heartdf$HeartDisease <- as.factor(heartdf$HeartDisease)

sapply(heartdf, class)

##            Age            Sex  ChestPainType      RestingBP    Cholesterol 
##      "integer"       "factor"       "factor"      "integer"      "integer" 
##      FastingBS     RestingECG          MaxHR ExerciseAngina        Oldpeak 
##       "factor"       "factor"      "integer"       "factor"      "numeric" 
##       ST_Slope   HeartDisease 
##       "factor"       "factor"

EXPLANATION:

Ada beberapa atribut data yang memiliki tipe data character sebenarnya bersifat kategorik, selanjutnya saya mengubahnya menjadi factor. Lalu, saya mengecek kembali perubahan tipe data pada dataset yang saya gunakan, tipe data sudah berubah.

colSums(is.na(heartdf))

##            Age            Sex  ChestPainType      RestingBP    Cholesterol 
##              0              0              0              0              0 
##      FastingBS     RestingECG          MaxHR ExerciseAngina        Oldpeak 
##              0              0              0              0              0 
##       ST_Slope   HeartDisease 
##              0              0

EXPLANATION:

Saya mencari nilai missing value pada dataset yang saya gunakan. Dapat dilihat, tidak ada missing value pada dataset.

sapply(heartdf, function(x) length(unique(x)))

##            Age            Sex  ChestPainType      RestingBP    Cholesterol 
##             50              2              4             67            222 
##      FastingBS     RestingECG          MaxHR ExerciseAngina        Oldpeak 
##              2              3            119              2             53 
##       ST_Slope   HeartDisease 
##              3              2

EXPLANATION:

• Saya mengecek data yang unik di tiap atribut pada dataset. Pada tipe data numeric, terdapat banyak data yang unik karena nilai angka di dalam nya berbeda - beda. Contohnya: variabel Age menunjukkan umur para pengidap penyakit jantung, tentunya tiap orang memiliki umur yang berbeda saat mengalami penyakit jantung.

• Sedangkan, pada tipe data factor hanya terdapat beberapa data yang unik, karena mewakilkan kategori yang ada dalam tiap atribut data. Contohnya: variabel HeartDisease, 2 kategori dalam HeartDisease menunjukkan seseorang mengalami penyakit jantung dan tidak sakit jantung.

ThreeSigma <- function(x, t = 3){
 mu <- mean(x, na.rm = TRUE)
 sig <- sd(x, na.rm = TRUE)
 if (sig == 0){
 message("All non-missing x-values are identical")
}
 up <- mu + t * sig
 down <- mu - t * sig
 out <- list(up = up, down = down)
 
 return(out)
 }

Hampel <- function(x, t = 3){
 mu <- median(x, na.rm = TRUE)
 sig <- mad(x, na.rm = TRUE)
 if (sig == 0){
 message("Hampel identifer implosion: MAD scale estimate is zero")
 }
 up <- mu + t * sig
 down <- mu - t * sig
 out <- list(up = up, down = down)
 
 return(out)
 }
   
BoxplotRule<- function(x, t = 1.5){
 xL <- quantile(x, na.rm = TRUE, probs = 0.25, names = FALSE)
 xU <- quantile(x, na.rm = TRUE, probs = 0.75, names = FALSE)
 Q <- xU - xL
 if (Q == 0){
 message("Boxplot rule implosion: interquartile distance is zero")
 }
 up <- xU + t * Q
 down <- xU - t * Q
 out <- list(up = up, down = down)
 
 return(out)
}   

ExtractDetails <- function(x, down, up){
 outClass <- rep("N", length(x))
 indexLo <- which(x < down)
 indexHi <- which(x > up)
 outClass[indexLo] <- "L"
 outClass[indexHi] <- "U"
 index <- union(indexLo, indexHi)
 values <- x[index]
 outClass <- outClass[index]
 nOut <- length(index)
 maxNom <- max(x[which(x <= up)])
 minNom <- min(x[which(x >= down)])
 outList <- list(nOut = nOut, lowLim = down,
 upLim = up, minNom = minNom,
 maxNom = maxNom, index = index,
 values = values,
 outClass = outClass)
 
 return(outList)
 }

FindOutliers <- function(x, t3 = 3, tH = 3, tb = 1.5){
 threeLims <- ThreeSigma(x, t = t3)
 HampLims <- Hampel(x, t = tH)
 boxLims <- BoxplotRule(x, t = tb)

 n <- length(x)
 nMiss <- length(which(is.na(x)))

 threeList <- ExtractDetails(x, threeLims$down, threeLims$up)
 HampList <- ExtractDetails(x, HampLims$down, HampLims$up)
 boxList <- ExtractDetails(x, boxLims$down, boxLims$up)

 sumFrame <- data.frame(method = "ThreeSigma", n = n,
 nMiss = nMiss, nOut = threeList$nOut,
 lowLim = threeList$lowLim,
 upLim = threeList$upLim,
 minNom = threeList$minNom,
 maxNom = threeList$maxNom)
 upFrame <- data.frame(method = "Hampel", n = n,
 nMiss = nMiss, nOut = HampList$nOut,
 lowLim = HampList$lowLim,
 upLim = HampList$upLim,
 minNom = HampList$minNom,
 maxNom = HampList$maxNom)
 sumFrame <- rbind.data.frame(sumFrame, upFrame)
 upFrame <- data.frame(method = "BoxplotRule", n = n,
 nMiss = nMiss, nOut = boxList$nOut,
 lowLim = boxList$lowLim,
 upLim = boxList$upLim,
 minNom = boxList$minNom,
 maxNom = boxList$maxNom)
 sumFrame <- rbind.data.frame(sumFrame, upFrame)

 threeFrame <- data.frame(index = threeList$index,
 values = threeList$values,
 type = threeList$outClass)
 HampFrame <- data.frame(index = HampList$index,
 values = HampList$values,
 type = HampList$outClass)
 boxFrame <- data.frame(index = boxList$index,
 values = boxList$values,
 type = boxList$outClass)
 outList <- list(summary = sumFrame, threeSigma = threeFrame,
 Hampel = HampFrame, boxplotRule = boxFrame)
 
 return(outList)
}

outlier_oldpeak <- FindOutliers(heartdf$Oldpeak)
outlier_oldpeak$summary

##        method   n nMiss nOut    lowLim    upLim minNom maxNom
## 1  ThreeSigma 918     0    7 -2.312347 4.087074   -2.0    4.0
## 2      Hampel 918     0   26 -2.068680 3.268680   -2.0    3.2
## 3 BoxplotRule 918     0   23 -0.750000 3.750000   -0.7    3.7

outlier_restBP <- FindOutliers(heartdf$RestingBP)
outlier_restBP$summary

##        method   n nMiss nOut    lowLim   upLim minNom maxNom
## 1  ThreeSigma 918     0    8  76.85405 187.939     80    185
## 2      Hampel 918     0   25  85.52200 174.478     92    174
## 3 BoxplotRule 918     0   80 110.00000 170.000    110    170

outlier_maxHR <- FindOutliers(heartdf$MaxHR)
outlier_maxHR$summary

##        method   n nMiss nOut    lowLim    upLim minNom maxNom
## 1  ThreeSigma 918     0    1  60.42837 213.1904     63    202
## 2      Hampel 918     0    0  57.93960 218.0604     60    202
## 3 BoxplotRule 918     0   85 102.00000 210.0000    102    202

outlier_cholesterol <- FindOutliers(heartdf$Cholesterol)
outlier_cholesterol$summary

##        method   n nMiss nOut    lowLim    upLim minNom maxNom
## 1  ThreeSigma 918     0    3 -129.3529 526.9520      0    518
## 2      Hampel 918     0  180   18.4012 427.5988     85    417
## 3 BoxplotRule 918     0  192  126.3750 407.6250    129    407

outlier_Age <- FindOutliers(heartdf$Age)
outlier_Age$summary

##        method   n nMiss nOut   lowLim    upLim minNom maxNom
## 1  ThreeSigma 918     0    0 25.21304 81.80874     28     77
## 2      Hampel 918     0    0 22.86540 85.13460     28     77
## 3 BoxplotRule 918     0   93 40.50000 79.50000     41     77

EXPLANATION:

• Saya mencari data outlier pada tiap variabel numeric menggunakan metode three sigma, hampel identifier, dan boxplot rule. Dari ketiga metode, saya memilih boxplot rule untuk menangani nilai outlier yang ditemukan pada tiap variabel karena menghasilkan data outlier paling banyak dari 2 metode lainnya.

• Dalam logistic regression, adanya outlier tidak terlalu berpengaruh terhadap model regresi. Tetapi, data outlier yang ekstrim dapat berpengaruh terhadap model regresi yang akan dibuat. Sehingga, banyak data outlier yang bernilai sekitar 10% dari total data tidak akan saya gunakan.

NoOutlierDF <- subset(heartdf, heartdf$RestingBP < 170.000 & heartdf$RestingBP > 110.000)  
NoOutlierDF <- subset(heartdf, heartdf$MaxHR < 210.0000 & heartdf$MaxHR > 102.0000)
NoOutlierDF <- subset(heartdf, heartdf$Cholesterol < 407.6250 & heartdf$Cholesterol > 126.3750)
NoOutlierDF <- subset(heartdf, heartdf$Age < 79.5000 & heartdf$Age > 40.5000)

dim(NoOutlierDF)

## [1] 825  12

EXPLANATION:

Saya mengambil data yang tidak mengandung outlier menggunakan function subset(). Range data yang saya gunakan menggunakan batas atas dan bawah dari perhitungan boxplot rule. Selanjutnya saya cek dimensi dari dataframe yang tidak mengandung outlier, data outlier sudah tidak termasuk dalam dataframe.

Univariate

heartdf_numeric <- (NoOutlierDF[sapply(NoOutlierDF, is.numeric)])
hist.data.frame(heartdf_numeric)

EXPLANATION:

Saya melalukan visualisasi terhadap variabel yang bersifat numeric. Dari visualisasi diatas dapat disimpulkan :

1. Variabel Age (umur) tidak terdistribusi secara normal, visualisasinya tidak simetris.
2. Variabel RestingBP (tekanan darah istirahat) memiliki beberapa data yang rentangnya sangat jauh dari data lainnya, sehingga tidak terdistribusi secara normal.
3. Variabel Cholesterol (kolesterol) hampir terdistribusi secara normal, visualisasi mengarah ke bentuk bell curve walaupun tidak sempurna. Selain itu, masih ada beberapa pencilan data pada variabel Cholesterol membuat tidak sepenuhnya terdistribusi secara normal.
4. Variabel MaxHR (detak jantung maksimum) hampir terdistribusi secara normal, visualisasi hampir berbentuk bell curve tetapi adanya beberapa yang tidak simetris membuat visualisasi nya tidak sepenuhnya baik.
5. Variabel Oldpeak tidak terdistribusi secara normal, visualisasi tidak simetris. Selain itu, ada beberapa pencilan data pada variabel nya.

heartdf_factor <- (NoOutlierDF[sapply(NoOutlierDF, is.factor)])
hist.data.frame(heartdf_factor)

EXPLANATION:

Saya melalukan visualisasi terhadap variabel yang bersifat categorical. Dari visualisasi diatas dapat disimpulkan :

1. Variabel Sex di dominasi oleh “Male” (pria) sebanyak lebih dari 600 orang. Sebagian lain terdiri dari “Female” (perempuan) kurang lebih 100 orang.
2. Variabel ChestPainType memiliki data terbanyak pada tipe nyeri dada “ASY” (lebih dari 400 orang). Selanjutnya, ada tipe nyeri dada “NAP” & “ATA” (sekitar 100 orang), serta tipe nyeri dada TA (kurang dari 100 orang) dengan jumlah paling sedikit.
3. Variabel FastingBS yaitu gula darah puasa (1 = tinggi, 0 = normal atau rendah) memiliki 2 kategori. Kategori tinggi (sebanyak 200 orang) dan selain tinggi (gula darah normal atau rendah) sebanyak kurang lebih 600 orang.
4. Variabel RestingECG yaitu elektrokardiogram istirahat dengan kategori ST (kelainan gelombang ST-T) & LVH (kemungkinan hipertrofi ventrikel kiri) kurang lebih 200 orang. Sedangkan data terbanyak dengan kategori normal sebanyak lebih dari 450 orang.
5. Variabel ExerciseAngina, yaitu penyakit angina akibat jantung beraktivitas lebih keras (seperti olahraga) dialami oleh kurang dari 360 orang. Sebagian besar lain tidak mengalami angina, lebih dari 460 orang.
6. Variabel ST_Slope, yaitu puncak kemiringan saat latihan ST, dengan kategori down (menurun) sebanyak kurang dari 100 orang, up (menanjak) sebanyak kurang lebih 300 orang, dan flat (datar) sebanyak lebih dari 400 orang.
7. Variabel HeartDisease, yaitu sakit jantung memiliki 2 kategori (1 = penderita sakit jantung, 0 = tidak sakit jantung). Orang yang mengalami penyakit jantung (sekitar 480 orang) lebih banyak dari yang tidak sakit jantung. Sedangkan yang tidak mengalami penyakit jantung kurang dari 360 orang.

Bivariate

ggplot(NoOutlierDF, aes(x = HeartDisease, fill = Sex)) + 
  geom_bar(position = "fill") + 
  ylab("proportion")

EXPLANATION:

Dari visualisasi yang saya lakukan, dapat disimpulkan penyakit jantung lebih banyak dialami oleh laki - laki.

ggplot(NoOutlierDF, aes(x = HeartDisease, fill = ChestPainType)) +
  geom_bar(position = "dodge")

EXPLANATION:

Penderita penyakit jantung paling sering mengalami nyeri dada dengan tipe “ASY”. Sedangkan orang yang tidak mengalami sakit jantung paling banyak mengalami nyeri dada dengan tipe “ATA”.

ggplot(NoOutlierDF, aes(x = HeartDisease, fill = FastingBS)) + 
  geom_bar(position = "fill") + 
  ylab("proportion")

EXPLANATION:

Penderita penyakit jantung cenderung tidak memiliki gula darah puasa yang tinggi (1: gula darah puasa tinggi, 0: normal atau rendah). Gula darah puasa yang tinggi tidak terlalu berpengaruh terhadap penderita penyakit jantung.

ggplot(NoOutlierDF, aes(x = HeartDisease, fill = RestingECG)) + 
  geom_bar(position = "fill") + 
  ylab("proportion")

EXPLANATION:

Penderita penyakit jantung kebanyakan memiliki elektrokardiogram istirahat yang normal.

ggplot(NoOutlierDF, aes(x = HeartDisease, fill = ExerciseAngina)) + 
  geom_bar(position = "fill") + 
  ylab("proportion")

EXPLANATION:

Penderita penyakit jantung kebanyakan mengalami sakit angina yang disebabkan saat jantung sedang bekerja lebih keras. Sedangkan orang yang tidak mengalami penyakit jantung kebanyakan tidak mengalami sakit angina (nyeri di dada).

ggplot(NoOutlierDF, aes(x = HeartDisease, fill = ST_Slope)) + 
  geom_bar(position = "fill") + 
  ylab("proportion")

EXPLANATION:

Penderita penyakit jantung sebagian besar memiliki kemiringan yang datar (Flat) saat latihan ST.

Multivariate

heartdf_numeric <- (NoOutlierDF[sapply(NoOutlierDF, is.numeric)])
rcorr(as.matrix(heartdf_numeric), type = "spearman")

##               Age RestingBP Cholesterol MaxHR Oldpeak
## Age          1.00      0.24       -0.06 -0.29    0.25
## RestingBP    0.24      1.00        0.11 -0.09    0.18
## Cholesterol -0.06      0.11        1.00  0.22    0.06
## MaxHR       -0.29     -0.09        0.22  1.00   -0.16
## Oldpeak      0.25      0.18        0.06 -0.16    1.00
## 
## n= 825 
## 
## 
## P
##             Age    RestingBP Cholesterol MaxHR  Oldpeak
## Age                0.0000    0.0705      0.0000 0.0000 
## RestingBP   0.0000           0.0019      0.0121 0.0000 
## Cholesterol 0.0705 0.0019                0.0000 0.1007 
## MaxHR       0.0000 0.0121    0.0000             0.0000 
## Oldpeak     0.0000 0.0000    0.1007      0.0000

EXPLANATION:

Untuk melihat berbagai keterkaitan antar variabel numeric saya menggunakan function rcorr(). Dapat dilihat hubungan terkuat ada pada variabel Oldpeak dengan Age.

1B. Data Preparation

DFTrain <- subset(NoOutlierDF, select = c(1, 2, 3, 5, 7, 8, 9, 10, 12))

EXPLANATION:

Saya memilih kolom dari variabel yang akan saya gunakan sebagai training data. Saya tidak menggunakan variabel RestingBP (tekanan darah istirahat), FastingBS (gula darah puasa), dan ST_Slope (kemiringan saat latihan ST) karena pengaruhnya cukup rendah terhadap HeartDisease. Hal ini dapat dilihat dari korelasi variabelnya dengan HeartDisease.

dim(DFTrain)

## [1] 825   9

head(DFTrain)

##    Age Sex ChestPainType Cholesterol RestingECG MaxHR ExerciseAngina Oldpeak
## 2   49   F           NAP         180     Normal   156              N     1.0
## 4   48   F           ASY         214     Normal   108              Y     1.5
## 5   54   M           NAP         195     Normal   122              N     0.0
## 7   45   F           ATA         237     Normal   170              N     0.0
## 8   54   M           ATA         208     Normal   142              N     0.0
## 10  48   F           ATA         284     Normal   120              N     0.0
##    HeartDisease
## 2             1
## 4             1
## 5             0
## 7             0
## 8             0
## 10            0

colSums(is.na(DFTrain))

##            Age            Sex  ChestPainType    Cholesterol     RestingECG 
##              0              0              0              0              0 
##          MaxHR ExerciseAngina        Oldpeak   HeartDisease 
##              0              0              0              0

EXPLANATION:

Saya mengecek dataframe dari training set yang saya buat. Terdapat 825 baris dan 9 atribut dengan tipe data factor dan numeric. Tidak ada missing value pada training set.

Split Data Into Training and Validation Set

0.8 * 825

## [1] 660

TrainingSet <- DFTrain[1:660,]
ValidationSet <- DFTrain[661:825,]

EXPLANATION:

Rasio pembagian data training dan validation set yang saya lakukan:

- 80% data untuk training set

- 30% data untuk validation set

1C. Modelling

Saya ingin membuat logistic regression model dengan menggunakan HeartDisease sebagai target variable.

logistic_Fit1 <- glm(HeartDisease~.,family = binomial(link = "logit"), TrainingSet)
summary(logistic_Fit1)

## 
## Call:
## glm(formula = HeartDisease ~ ., family = binomial(link = "logit"), 
##     data = TrainingSet)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6725  -0.5522   0.3097   0.5734   2.6125  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       0.796353   1.368396   0.582 0.560593    
## Age               0.023297   0.015523   1.501 0.133411    
## SexM              1.065015   0.293760   3.625 0.000288 ***
## ChestPainTypeATA -2.003477   0.334382  -5.992 2.08e-09 ***
## ChestPainTypeNAP -1.010020   0.264324  -3.821 0.000133 ***
## ChestPainTypeTA  -0.850393   0.479968  -1.772 0.076433 .  
## Cholesterol      -0.004148   0.001028  -4.033 5.51e-05 ***
## RestingECGNormal -0.189821   0.312292  -0.608 0.543299    
## RestingECGST     -0.193754   0.355704  -0.545 0.585955    
## MaxHR            -0.014474   0.005085  -2.846 0.004421 ** 
## ExerciseAnginaY   1.311811   0.252232   5.201 1.98e-07 ***
## Oldpeak           0.543238   0.125121   4.342 1.41e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 886.73  on 659  degrees of freedom
## Residual deviance: 539.50  on 648  degrees of freedom
## AIC: 563.5
## 
## Number of Fisher Scoring iterations: 5

logistic_Fit2 <- glm(HeartDisease~ Sex + Cholesterol + ChestPainType + ExerciseAngina + Oldpeak, family = binomial(link = "logit"), TrainingSet)
summary(logistic_Fit2)

## 
## Call:
## glm(formula = HeartDisease ~ Sex + Cholesterol + ChestPainType + 
##     ExerciseAngina + Oldpeak, family = binomial(link = "logit"), 
##     data = TrainingSet)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7524  -0.5279   0.3176   0.6110   2.6892  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       0.098723   0.367622   0.269 0.788280    
## SexM              1.092954   0.289430   3.776 0.000159 ***
## Cholesterol      -0.004819   0.000988  -4.877 1.07e-06 ***
## ChestPainTypeATA -2.150167   0.324994  -6.616 3.69e-11 ***
## ChestPainTypeNAP -1.076439   0.258909  -4.158 3.22e-05 ***
## ChestPainTypeTA  -0.825548   0.466356  -1.770 0.076692 .  
## ExerciseAnginaY   1.470573   0.244385   6.017 1.77e-09 ***
## Oldpeak           0.551388   0.122322   4.508 6.55e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 886.73  on 659  degrees of freedom
## Residual deviance: 554.79  on 652  degrees of freedom
## AIC: 570.79
## 
## Number of Fisher Scoring iterations: 5

EXPLANATION:

Saya membuat 2 model, pertama logistic_Fit1 dengan menggunakan semua atribut. Kemudian, saya membuat model logistic_Fit2 dengan beberapa variabel yang memiliki P-value paling rendah dibandingkan variabel lainnya.

head(ValidationSet)

##     Age Sex ChestPainType Cholesterol RestingECG MaxHR ExerciseAngina Oldpeak
## 741  54   F           NAP         201     Normal   163              N     0.0
## 742  62   M           ASY         267     Normal    99              Y     1.8
## 743  52   F           NAP         196        LVH   169              N     0.1
## 744  52   M           ATA         201     Normal   158              N     0.8
## 745  60   M           ASY         230     Normal   160              Y     1.4
## 746  63   F           ASY         269     Normal   169              Y     1.8
##     HeartDisease
## 741            0
## 742            1
## 743            0
## 744            0
## 745            1
## 746            1

predict_Fit1 <- predict(logistic_Fit1, newdata = subset(ValidationSet, select = c(1, 2, 3, 4, 5, 6, 7, 8, 9)), type = "response")
predict_Fit2 <- predict(logistic_Fit2, newdata = subset(ValidationSet, select = c(1, 2, 3, 4, 5, 6,7, 8, 9)), type = "response")

EXPLANATION:

Saya ingin menguji model yang saya buat, yaitu logistic_Fit1 dan logistic_Fit2 terhadap validation set. Kemudian saya memilih kolom pada validation set sesuai model yang saya buat di dalam variabel predict_Fit1.

evaluation_Fit1 <-  prediction(predict_Fit1, ValidationSet$HeartDisease)
prf_Fit1 <- performance(evaluation_Fit1, measure = "tpr", x.measure = "fpr")
plot(prf_Fit1)

evaluation_Fit2 <-  prediction(predict_Fit2, ValidationSet$HeartDisease)
prf_Fit2 <- performance(evaluation_Fit2, measure = "tpr", x.measure = "fpr")
plot(prf_Fit2)

EXPLANATION:

Kurva ROC pada kedua model (Fit1 & Fit2) menunjukkan hasil yang baik. Lengkungan kurva terakhir (paling kanan) menunjukkan nilai diatas 0.8. Nilai tersebut tergolong baik, karena semakin mendekati 1.

Check AUC Score

# check AUC score validation set
auc1 <- performance(evaluation_Fit1, measure = "auc")
auc1 <- auc1@y.values[[1]]
auc1

## [1] 0.8872679

# check AUC score training set
predictFit1_train <- predict(logistic_Fit1, newdata = subset(TrainingSet, select = c(1, 2, 3, 4, 5, 6, 7, 8, 9)) ,type = "response")
evaluation_Fit1 <-  prediction(predictFit1_train, TrainingSet$HeartDisease)

auc1train <- performance(evaluation_Fit1, measure = "auc")
auc1train <- auc1train@y.values[[1]]
auc1train

## [1] 0.887232

EXPLANATION:

Nilai AUC Score (Fit1) pada validation dan training set tidak terlalu jauh. Hal ini menandakan model yang dibuat tidak bias dan overfit.

auc2 <- performance(evaluation_Fit2, measure = "auc")
auc2 <- auc2@y.values[[1]]
auc2

## [1] 0.8763631

EXPLANATION:

Nilai AUC Score pada Fit2 cukup baik, tetapi saya menggunakan model dengan nilai AUC score yang paling tinggi, yaitu Fit 1.

2. Decision Tree

training_DTree <- DFTrain[1:660,]
validation_DTree <- DFTrain[661:825,]

EXPLANATION:

Untuk decision tree, saya juga membagi data dengan perbandingan:

- 80% data untuk training set

- 20% data untuk validation set

modelDTree <- rpart(HeartDisease~., data = training_DTree)
modelDTree

## n= 660 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 660 262 1 (0.39696970 0.60303030)  
##    2) ChestPainType=ATA,NAP,TA 274  90 0 (0.67153285 0.32846715)  
##      4) ExerciseAngina=N 214  47 0 (0.78037383 0.21962617)  
##        8) Cholesterol>=119.5 186  29 0 (0.84408602 0.15591398) *
##        9) Cholesterol< 119.5 28  10 1 (0.35714286 0.64285714) *
##      5) ExerciseAngina=Y 60  17 1 (0.28333333 0.71666667)  
##       10) MaxHR>=150.5 7   0 0 (1.00000000 0.00000000) *
##       11) MaxHR< 150.5 53  10 1 (0.18867925 0.81132075) *
##    3) ChestPainType=ASY 386  78 1 (0.20207254 0.79792746)  
##      6) Oldpeak< 0.65 144  52 1 (0.36111111 0.63888889)  
##       12) Cholesterol>=42.5 93  44 0 (0.52688172 0.47311828)  
##         24) ExerciseAngina=N 64  24 0 (0.62500000 0.37500000) *
##         25) ExerciseAngina=Y 29   9 1 (0.31034483 0.68965517) *
##       13) Cholesterol< 42.5 51   3 1 (0.05882353 0.94117647) *
##      7) Oldpeak>=0.65 242  26 1 (0.10743802 0.89256198) *

rpart.plot(modelDTree)

EXPLANATION:

Decision Tree dapat digunakan untuk memprediksi data baru yang belum memiliki label (keputusan akhir) caranya dengan menyesuaikan kriteria data sesuai branch yang ada pada Decision Tree. (Hasil akhir 1: memprediksi mengalami penyakit jantung, 0: memprediksi orang yang tidak mengalami penyakit jantung).

Important Variable

modelDTree$variable.importance

##  ChestPainType ExerciseAngina        Oldpeak          MaxHR    Cholesterol 
##      70.635282      47.524852      29.576620      26.705072      25.972802 
##            Sex            Age 
##      10.827306       8.357772

EXPLANATION:

Pada decision tree yang sudah dibuat, variabel ChestPainType (root) merupakan variable terpenting dalam tree. Karena nilainya paling besar dari variabel lainnya.

predictionDTree <- predict(modelDTree, newdata = subset(validation_DTree, select = c(1, 2, 3, 4, 5, 6,7, 8, 9)), type = "class")

Conclusion

table(predictionDTree, validation_DTree$HeartDisease)

##                
## predictionDTree  0  1
##               0 77 24
##               1 10 54

EXPLANATION:

Tabel prediction diatas menunjukkan hasil:

1. true-negative sebanyak 77 (ada 77 orang yang diprediksi tidak sakit jantung dan benar tidak sakit jantung)
2. true-positive sebanyak 54 (ada 54 orang yang diprediksi sakit jantung dan benar sakit jantung)
3. false-negative sebanyak 10 (ada 10 orang yang diprediksi tidak sakit jantung tetapi sebenarnya sakit jantung / prediksi salah)
4. false-positive sebanyak 24 (ada 24 orang yang diprediksi sakit jantung tetapi sebenarnya tidak sakit jantung / prediksi salah)

correct_clasif <- table(predictionDTree, validation_DTree$HeartDisease)
sum(diag(correct_clasif))  #correct classification data

## [1] 131

EXPLANATION:

Ada 131 data yang hasil prediksinya benar.

sum(correct_clasif) - sum(diag(correct_clasif)) #incorrecct classification data

## [1] 34

EXPLANATION:

Ada 34 data yang hasil prediksinya salah.

Evaluation

# accuracy per category 
accuracy <- diag(correct_clasif) / rowSums(correct_clasif) * 100
accuracy

##        0        1 
## 76.23762 84.37500

EXPLANATION:

Nilai akurasi untuk memprediksi kategori penderita penyakit jantung sebesar 84%. Sedangkan, nilai akurasi untuk memprediksi orang yang tidak mengalami sakit jantung sebesar 76%.

overallAccuracy <- sum(diag(correct_clasif)) / sum(correct_clasif) * 100
overallAccuracy

## [1] 79.39394

EXPLANATION:

Nilai akurasi keseluruhan decision tree adalah 79%. Menurut saya, nilai tersebut cukup baik karena akurasi cukup mendekati angka 100%.