Heart Disease Risk Prediction
Jasmine Mutia Alifa
2022-07-11
Link Video Penjelasan : https://youtu.be/fKa1mHW8hNs
Load Data
library(readr)
library(ggplot2)
library(Hmisc)##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(dplyr)##
## 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, unionlibrary(ROCR)
library(rpart)
library(rpart.plot)
library(caret)## Loading required package: latticepath <- "D:/R/heartData.csv"
heartdf <- read.csv(path)1A. Exploratory Data Analysis
dim(heartdf)## [1] 918 12str(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 0EXPLANATION:
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 2EXPLANATION:
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.7outlier_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 170outlier_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 202outlier_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 407outlier_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 77EXPLANATION:
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 12EXPLANATION:
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 :
- Variabel
Age(umur) tidak terdistribusi secara normal, visualisasinya tidak simetris. - Variabel
RestingBP(tekanan darah istirahat) memiliki beberapa data yang rentangnya sangat jauh dari data lainnya, sehingga tidak terdistribusi secara normal. - Variabel
Cholesterol(kolesterol) hampir terdistribusi secara normal, visualisasi mengarah ke bentuk bell curve walaupun tidak sempurna. Selain itu, masih ada beberapa pencilan data pada variabelCholesterolmembuat tidak sepenuhnya terdistribusi secara normal. - 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. - Variabel
Oldpeaktidak 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 :
- Variabel
Sexdi dominasi oleh “Male” (pria) sebanyak lebih dari 600 orang. Sebagian lain terdiri dari “Female” (perempuan) kurang lebih 100 orang. - Variabel
ChestPainTypememiliki 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. - Variabel
FastingBSyaitu 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. - Variabel
RestingECGyaitu 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. - 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. - 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. - 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.0000EXPLANATION:
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 9head(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 0colSums(is.na(DFTrain))## Age Sex ChestPainType Cholesterol RestingECG
## 0 0 0 0 0
## MaxHR ExerciseAngina Oldpeak HeartDisease
## 0 0 0 0EXPLANATION:
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] 660TrainingSet <- 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)
##
## 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: 5logistic_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)
##
## 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: 5EXPLANATION:
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.
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
## 747 66 M ASY 212 LVH 132 Y 0.1
## 748 42 M ASY 226 Normal 178 N 0.0
## 749 64 M ASY 246 LVH 96 Y 2.2
## 750 54 M NAP 232 LVH 165 N 1.6
## 751 46 F NAP 177 LVH 160 Y 1.4
## 752 67 F NAP 277 Normal 172 N 0.0
## 753 56 M ASY 249 LVH 144 Y 1.2
## 755 57 M ASY 207 Normal 168 Y 0.0
## 756 64 M ASY 212 LVH 132 N 2.0
## 757 59 M ASY 271 LVH 182 N 0.0
## 758 50 M NAP 233 Normal 163 N 0.6
## 759 51 M TA 213 LVH 125 Y 1.4
## 760 54 M ATA 283 LVH 195 N 0.0
## 761 53 M ASY 282 Normal 95 Y 2.0
## 762 52 M ASY 230 Normal 160 N 0.0
## 764 58 M NAP 224 LVH 173 N 3.2
## 765 41 F NAP 268 LVH 172 Y 0.0
## 766 41 M NAP 250 Normal 179 N 0.0
## 767 50 F NAP 219 Normal 158 N 1.6
## 768 54 F NAP 267 LVH 167 N 0.0
## 769 64 F ASY 303 Normal 122 N 2.0
## 770 51 F NAP 256 LVH 149 N 0.5
## 771 46 F ATA 204 Normal 172 N 0.0
## 772 55 M ASY 217 Normal 111 Y 5.6
## 773 45 M ATA 308 LVH 170 N 0.0
## 774 56 M TA 193 LVH 162 N 1.9
## 775 66 F ASY 228 Normal 165 Y 1.0
## 777 62 F ASY 244 Normal 154 Y 1.4
## 778 55 M ATA 262 Normal 155 N 0.0
## 779 58 M ASY 259 LVH 130 Y 3.0
## 780 43 M ASY 211 Normal 161 N 0.0
## 781 64 F ASY 325 Normal 154 Y 0.0
## 782 50 F ASY 254 LVH 159 N 0.0
## 783 53 M NAP 197 LVH 152 N 1.2
## 784 45 F ASY 236 LVH 152 Y 0.2
## 785 65 M TA 282 LVH 174 N 1.4
## 786 69 M TA 234 LVH 131 N 0.1
## 787 69 M NAP 254 LVH 146 N 2.0
## 788 67 M ASY 299 LVH 125 Y 0.9
## 789 68 F NAP 211 LVH 115 N 1.5
## 791 62 F ASY 294 Normal 106 N 1.9
## 792 51 M ASY 298 Normal 122 Y 4.2
## 793 46 M NAP 231 Normal 147 N 3.6
## 794 67 M ASY 254 Normal 163 N 0.2
## 795 50 M NAP 196 Normal 163 N 0.0
## 796 42 M NAP 240 Normal 194 N 0.8
## 797 56 F ASY 409 LVH 150 Y 1.9
## 798 41 M ASY 172 LVH 158 N 0.0
## 799 42 F ASY 265 LVH 122 N 0.6
## 800 53 M NAP 246 LVH 173 N 0.0
## 801 43 M NAP 315 Normal 162 N 1.9
## 802 56 M ASY 184 LVH 105 Y 2.1
## 803 52 M ASY 233 Normal 147 N 0.1
## 804 62 F ASY 394 LVH 157 N 1.2
## 805 70 M NAP 269 Normal 112 Y 2.9
## 806 54 M ASY 239 Normal 160 N 1.2
## 807 70 M ASY 174 Normal 125 Y 2.6
## 808 54 M ATA 309 Normal 156 N 0.0
## 810 48 M NAP 255 Normal 175 N 0.0
## 811 55 F ATA 250 LVH 161 N 1.4
## 812 58 F ASY 248 LVH 122 N 1.0
## 813 54 F NAP 214 Normal 158 N 1.6
## 814 69 F TA 239 Normal 151 N 1.8
## 815 77 M ASY 304 LVH 162 Y 0.0
## 816 68 M NAP 277 Normal 151 N 1.0
## 817 58 M ASY 300 LVH 171 N 0.0
## 818 60 M ASY 258 LVH 141 Y 2.8
## 819 51 M ASY 299 Normal 173 Y 1.6
## 820 55 M ASY 289 LVH 145 Y 0.8
## 821 52 M TA 298 Normal 178 N 1.2
## 822 60 F NAP 318 Normal 160 N 0.0
## 823 58 M NAP 240 LVH 154 Y 0.6
## 824 64 M NAP 309 Normal 131 Y 1.8
## 826 59 M TA 288 LVH 159 N 0.2
## 827 51 M NAP 245 LVH 166 N 2.4
## 828 43 F NAP 213 Normal 165 N 0.2
## 829 58 M ASY 216 LVH 131 Y 2.2
## 831 41 F ATA 204 LVH 172 N 1.4
## 832 63 F NAP 252 LVH 172 N 0.0
## 833 51 M NAP 227 Normal 154 Y 0.0
## 834 54 M NAP 258 LVH 147 N 0.4
## 835 44 M ATA 220 Normal 170 N 0.0
## 836 54 M ASY 239 Normal 126 Y 2.8
## 837 65 M ASY 254 LVH 127 N 2.8
## 838 57 M NAP 168 Normal 174 N 1.6
## 839 63 M ASY 330 LVH 132 Y 1.8
## 841 41 M ATA 203 Normal 132 N 0.0
## 842 62 F NAP 263 Normal 97 N 1.2
## 843 43 F ASY 341 LVH 136 Y 3.0
## 844 58 F TA 283 LVH 162 N 1.0
## 845 52 M TA 186 LVH 190 N 0.0
## 846 61 F ASY 307 LVH 146 Y 1.0
## 848 45 M ASY 260 LVH 185 N 0.0
## 849 52 M ASY 255 Normal 161 Y 0.0
## 850 62 M NAP 231 Normal 146 N 1.8
## 851 62 F ASY 164 LVH 145 N 6.2
## 852 53 F ASY 234 LVH 160 N 0.0
## 853 43 M ASY 177 LVH 120 Y 2.5
## 854 47 M NAP 257 LVH 156 N 0.0
## 855 52 M ATA 325 Normal 172 N 0.2
## 856 68 M NAP 274 LVH 150 Y 1.6
## 858 53 F ASY 264 LVH 143 N 0.4
## 859 62 F ASY 268 LVH 160 N 3.6
## 860 51 F NAP 308 LVH 142 N 1.5
## 861 60 M ASY 253 Normal 144 Y 1.4
## 862 65 M ASY 248 LVH 158 N 0.6
## 863 65 F NAP 269 Normal 148 N 0.8
## 864 60 M NAP 185 LVH 155 N 3.0
## 865 60 M ASY 282 LVH 142 Y 2.8
## 866 54 M ASY 188 Normal 113 N 1.4
## 867 44 M ATA 219 LVH 188 N 0.0
## 868 44 M ASY 290 LVH 153 N 0.0
## 869 51 M NAP 175 Normal 123 N 0.6
## 870 59 M NAP 212 Normal 157 N 1.6
## 871 71 F ATA 302 Normal 162 N 0.4
## 872 61 M NAP 243 Normal 137 Y 1.0
## 873 55 M ASY 353 Normal 132 Y 1.2
## 874 64 M NAP 335 Normal 158 N 0.0
## 875 43 M ASY 247 Normal 171 N 1.5
## 876 58 F NAP 340 Normal 172 N 0.0
## 877 60 M ASY 206 LVH 132 Y 2.4
## 878 58 M ATA 284 LVH 160 N 1.8
## 879 49 M ATA 266 Normal 171 N 0.6
## 880 48 M ATA 229 Normal 168 N 1.0
## 881 52 M NAP 199 Normal 162 N 0.5
## 882 44 M ATA 263 Normal 173 N 0.0
## 883 56 F ATA 294 LVH 153 N 1.3
## 884 57 M ASY 192 Normal 148 N 0.4
## 885 67 M ASY 286 LVH 108 Y 1.5
## 886 53 F NAP 216 LVH 115 N 0.0
## 887 52 M NAP 223 Normal 169 N 0.0
## 888 43 M ASY 247 LVH 143 Y 0.1
## 889 52 M ASY 204 Normal 156 Y 1.0
## 890 59 M TA 204 Normal 162 N 0.8
## 891 64 M TA 227 LVH 155 N 0.6
## 892 66 F NAP 278 LVH 152 N 0.0
## 894 57 M ATA 232 LVH 164 N 0.0
## 895 58 F ASY 197 Normal 131 N 0.6
## 896 57 M ASY 335 Normal 143 Y 3.0
## 897 47 M NAP 253 Normal 179 N 0.0
## 898 55 F ASY 205 ST 130 Y 2.0
## 900 61 M ASY 203 Normal 161 N 0.0
## 901 58 M ASY 318 ST 140 N 4.4
## 902 58 F ASY 225 LVH 146 Y 2.8
## 903 58 M ATA 220 Normal 144 N 0.4
## 904 56 M ATA 221 LVH 163 N 0.0
## 905 56 M ATA 240 Normal 169 N 0.0
## 906 67 M NAP 212 LVH 150 N 0.8
## 907 55 F ATA 342 Normal 166 N 1.2
## 908 44 M ASY 169 Normal 144 Y 2.8
## 909 63 M ASY 187 LVH 144 Y 4.0
## 910 63 F ASY 197 Normal 136 Y 0.0
## 911 41 M ATA 157 Normal 182 N 0.0
## 912 59 M ASY 176 LVH 90 N 1.0
## 913 57 F ASY 241 Normal 123 Y 0.2
## 914 45 M TA 264 Normal 132 N 1.2
## 915 68 M ASY 193 Normal 141 N 3.4
## 916 57 M ASY 131 Normal 115 Y 1.2
## 917 57 F ATA 236 LVH 174 N 0.0
## HeartDisease
## 741 0
## 742 1
## 743 0
## 744 0
## 745 1
## 746 1
## 747 1
## 748 0
## 749 1
## 750 0
## 751 0
## 752 0
## 753 1
## 755 0
## 756 1
## 757 0
## 758 1
## 759 0
## 760 1
## 761 1
## 762 1
## 764 1
## 765 0
## 766 0
## 767 0
## 768 0
## 769 0
## 770 0
## 771 0
## 772 1
## 773 0
## 774 0
## 775 1
## 777 1
## 778 0
## 779 1
## 780 0
## 781 0
## 782 0
## 783 0
## 784 0
## 785 1
## 786 0
## 787 1
## 788 1
## 789 0
## 791 1
## 792 1
## 793 1
## 794 1
## 795 0
## 796 0
## 797 1
## 798 1
## 799 0
## 800 0
## 801 0
## 802 1
## 803 0
## 804 0
## 805 1
## 806 0
## 807 1
## 808 0
## 810 0
## 811 0
## 812 0
## 813 0
## 814 0
## 815 1
## 816 0
## 817 1
## 818 1
## 819 1
## 820 1
## 821 0
## 822 0
## 823 0
## 824 1
## 826 1
## 827 0
## 828 0
## 829 1
## 831 0
## 832 0
## 833 0
## 834 0
## 835 0
## 836 1
## 837 1
## 838 0
## 839 1
## 841 0
## 842 1
## 843 1
## 844 0
## 845 0
## 846 1
## 848 0
## 849 1
## 850 0
## 851 1
## 852 0
## 853 1
## 854 0
## 855 0
## 856 1
## 858 0
## 859 1
## 860 0
## 861 1
## 862 1
## 863 0
## 864 1
## 865 1
## 866 1
## 867 0
## 868 1
## 869 0
## 870 0
## 871 0
## 872 0
## 873 1
## 874 1
## 875 0
## 876 0
## 877 1
## 878 1
## 879 0
## 880 1
## 881 0
## 882 0
## 883 0
## 884 0
## 885 1
## 886 0
## 887 0
## 888 1
## 889 1
## 890 1
## 891 0
## 892 0
## 894 1
## 895 0
## 896 1
## 897 0
## 898 1
## 900 1
## 901 1
## 902 1
## 903 0
## 904 0
## 905 0
## 906 1
## 907 0
## 908 1
## 909 1
## 910 1
## 911 0
## 912 1
## 913 1
## 914 1
## 915 1
## 916 1
## 917 1predict_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 predictFit1.
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.887232EXPLANATION:
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.8763631EXPLANATION:
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.357772EXPLANATION:
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 54EXPLANATION:
Tabel prediction diatas menunjukkan hasil:
- true-negative sebanyak 77 (ada 77 orang yang diprediksi tidak sakit jantung dan benar tidak sakit jantung)
- true-positive sebanyak 54 (ada 54 orang yang diprediksi sakit jantung dan benar sakit jantung)
- false-negative sebanyak 10 (ada 10 orang yang diprediksi tidak sakit jantung tetapi sebenarnya sakit jantung / prediksi salah)
- 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] 131EXPLANATION:
Ada 131 data yang hasil prediksinya benar.
sum(correct_clasif) - sum(diag(correct_clasif)) #incorrecct classification data## [1] 34EXPLANATION:
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.37500EXPLANATION:
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.39394EXPLANATION:
Nilai akurasi keseluruhan decision tree adalah 79%. Menurut saya, nilai tersebut cukup baik karena akurasi cukup mendekati angka 100%.