2. Data Preparation
(top)
### column "AHD" is the respoinse variable.
# table of the respoonse variable
table(Heart$AHD)
##
## No Yes
## 164 139
print(Heart, width=1000) # show all columns
## # A tibble: 303 × 15
## index Age Sex ChestPain RestBP Chol Fbs RestECG MaxHR ExAng Oldpeak
## <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 63 1 typical 145 233 1 2 150 0 2.3
## 2 2 67 1 asymptomatic 160 286 0 2 108 1 1.5
## 3 3 67 1 asymptomatic 120 229 0 2 129 1 2.6
## 4 4 37 1 nonanginal 130 250 0 0 187 0 3.5
## 5 5 41 0 nontypical 130 204 0 2 172 0 1.4
## 6 6 56 1 nontypical 120 236 0 0 178 0 0.8
## 7 7 62 0 asymptomatic 140 268 0 2 160 0 3.6
## 8 8 57 0 asymptomatic 120 354 0 0 163 1 0.6
## 9 9 63 1 asymptomatic 130 254 0 2 147 0 1.4
## 10 10 53 1 asymptomatic 140 203 1 2 155 1 3.1
## Slope Ca Thal AHD
## <dbl> <dbl> <chr> <chr>
## 1 3 0 fixed No
## 2 2 3 normal Yes
## 3 2 2 reversable Yes
## 4 3 0 normal No
## 5 1 0 normal No
## 6 1 0 normal No
## 7 3 2 normal Yes
## 8 1 0 normal No
## 9 2 1 reversable Yes
## 10 3 0 reversable Yes
## # ℹ 293 more rows
# Rename medv column as resp.
Heart2 <- Heart %>%
select(-"index") %>% # remove col named "index"
rename(resp=AHD) %>% # rename the column
relocate(resp) # move "resp" to 1st column
print(Heart2, width=1000)
## # A tibble: 303 × 14
## resp Age Sex ChestPain RestBP Chol Fbs RestECG MaxHR ExAng Oldpeak
## <chr> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 No 63 1 typical 145 233 1 2 150 0 2.3
## 2 Yes 67 1 asymptomatic 160 286 0 2 108 1 1.5
## 3 Yes 67 1 asymptomatic 120 229 0 2 129 1 2.6
## 4 No 37 1 nonanginal 130 250 0 0 187 0 3.5
## 5 No 41 0 nontypical 130 204 0 2 172 0 1.4
## 6 No 56 1 nontypical 120 236 0 0 178 0 0.8
## 7 Yes 62 0 asymptomatic 140 268 0 2 160 0 3.6
## 8 No 57 0 asymptomatic 120 354 0 0 163 1 0.6
## 9 Yes 63 1 asymptomatic 130 254 0 2 147 0 1.4
## 10 Yes 53 1 asymptomatic 140 203 1 2 155 1 3.1
## Slope Ca Thal
## <dbl> <dbl> <chr>
## 1 3 0 fixed
## 2 2 3 normal
## 3 2 2 reversable
## 4 3 0 normal
## 5 1 0 normal
## 6 1 0 normal
## 7 3 2 normal
## 8 1 0 normal
## 9 2 1 reversable
## 10 3 0 reversable
## # ℹ 293 more rows
# Turn class of some columns into <factor>
Heart2 <- Heart2 %>%
mutate(resp = as.factor(resp), # Turn these columns to <factor> instead of <double>
Sex = as.factor(Sex),
ChestPain= as.factor(ChestPain),
Fbs = as.factor(Fbs),
RestECG = as.factor(RestECG),
ExAng = as.factor(ExAng),
Thal = as.factor(Thal))
print(Heart2, width=1000)
## # A tibble: 303 × 14
## resp Age Sex ChestPain RestBP Chol Fbs RestECG MaxHR ExAng Oldpeak
## <fct> <dbl> <fct> <fct> <dbl> <dbl> <fct> <fct> <dbl> <fct> <dbl>
## 1 No 63 1 typical 145 233 1 2 150 0 2.3
## 2 Yes 67 1 asymptomatic 160 286 0 2 108 1 1.5
## 3 Yes 67 1 asymptomatic 120 229 0 2 129 1 2.6
## 4 No 37 1 nonanginal 130 250 0 0 187 0 3.5
## 5 No 41 0 nontypical 130 204 0 2 172 0 1.4
## 6 No 56 1 nontypical 120 236 0 0 178 0 0.8
## 7 Yes 62 0 asymptomatic 140 268 0 2 160 0 3.6
## 8 No 57 0 asymptomatic 120 354 0 0 163 1 0.6
## 9 Yes 63 1 asymptomatic 130 254 0 2 147 0 1.4
## 10 Yes 53 1 asymptomatic 140 203 1 2 155 1 3.1
## Slope Ca Thal
## <dbl> <dbl> <fct>
## 1 3 0 fixed
## 2 2 3 normal
## 3 2 2 reversable
## 4 3 0 normal
## 5 1 0 normal
## 6 1 0 normal
## 7 3 2 normal
## 8 1 0 normal
## 9 2 1 reversable
## 10 3 0 reversable
## # ℹ 293 more rows
Heart2 %>% select(resp) %>% table
## resp
## No Yes
## 164 139
# note that level is "No/Yes", not "Yes/No". <---- Important!
3. Routine Exploratory Analysis (class of resp should be fct)
(top)
Orig <- Heart2
#- Check for N/A in data. Remove if there's any.
summary(Orig)
## resp Age Sex ChestPain RestBP
## No :164 Min. :29.00 0: 97 asymptomatic:144 Min. : 94.0
## Yes:139 1st Qu.:48.00 1:206 nonanginal : 86 1st Qu.:120.0
## Median :56.00 nontypical : 50 Median :130.0
## Mean :54.44 typical : 23 Mean :131.7
## 3rd Qu.:61.00 3rd Qu.:140.0
## Max. :77.00 Max. :200.0
##
## Chol Fbs RestECG MaxHR ExAng Oldpeak
## Min. :126.0 0:258 0:151 Min. : 71.0 0:204 Min. :0.00
## 1st Qu.:211.0 1: 45 1: 4 1st Qu.:133.5 1: 99 1st Qu.:0.00
## Median :241.0 2:148 Median :153.0 Median :0.80
## Mean :246.7 Mean :149.6 Mean :1.04
## 3rd Qu.:275.0 3rd Qu.:166.0 3rd Qu.:1.60
## Max. :564.0 Max. :202.0 Max. :6.20
##
## Slope Ca Thal
## Min. :1.000 Min. :0.0000 fixed : 18
## 1st Qu.:1.000 1st Qu.:0.0000 normal :166
## Median :2.000 Median :0.0000 reversable:117
## Mean :1.601 Mean :0.6722 NA's : 2
## 3rd Qu.:2.000 3rd Qu.:1.0000
## Max. :3.000 Max. :3.0000
## NA's :4
## [1] 303 14
## [1] 6
# If there is na in the data, run below
dim(Orig)
## [1] 303 14
Orig <- Orig %>% na.omit()
dim(Orig)
## [1] 297 14
##----------
##- Correlation Check
# Pick columns that are numeric
library(corrplot) # install.packages("corrplot")
Orig_num <- Orig %>% select_if(is.numeric)
cor(Orig_num)
## Age RestBP Chol MaxHR Oldpeak
## Age 1.0000000 0.29047626 2.026435e-01 -3.945629e-01 0.19712262
## RestBP 0.2904763 1.00000000 1.315357e-01 -4.910766e-02 0.19124314
## Chol 0.2026435 0.13153571 1.000000e+00 -7.456799e-05 0.03859579
## MaxHR -0.3945629 -0.04910766 -7.456799e-05 1.000000e+00 -0.34763997
## Oldpeak 0.1971226 0.19124314 3.859579e-02 -3.476400e-01 1.00000000
## Slope 0.1594047 0.12117205 -9.215240e-03 -3.893067e-01 0.57903735
## Ca 0.3622103 0.09795376 1.159446e-01 -2.687270e-01 0.29445228
## Slope Ca
## Age 0.15940474 0.36221034
## RestBP 0.12117205 0.09795376
## Chol -0.00921524 0.11594459
## MaxHR -0.38930674 -0.26872698
## Oldpeak 0.57903735 0.29445228
## Slope 1.00000000 0.10976112
## Ca 0.10976112 1.00000000
corrplot::corrplot(cor(Orig_num))
corrplot::corrplot(cor(Orig_num), method="number")
##----------
##- Visualization
# pairs(Orig) equivalent
library(GGally) # install.packages("GGally")
GGally::ggpairs(Orig[, c(2:7, 1)], aes(color=resp, alpha=1))
GGally::ggpairs(Orig[, c(8:14, 1)], aes(color=resp, alpha=1))
###----------
###- Chi-sq test of association for each column
# Pick columns that has more than 1 unique value
list_cols <- Orig %>% summarise_all(function(x) length(unique(x))) %>%
gather() %>% filter(value>1) %>% pull(key)
# Apply chisq.test to those columns
ChiSq.pval <- Orig %>% select(list_cols) %>%
summarise_all(funs(chisq.test(., Orig$resp)$p.value))
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(list_cols)
##
## # Now:
## data %>% select(all_of(list_cols))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
##
## # Simple named list: list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
##
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: There were 6 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `Age = chisq.test(Age, Orig$resp)$p.value`.
## Caused by warning in `chisq.test()`:
## ! Chi-squared approximation may be incorrect
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 5 remaining warnings.
ChiSq.pval <- Orig %>% summarise_all(funs(chisq.test(., Orig$resp)$p.value))
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
##
## # Simple named list: list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
##
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## There were 6 warnings in `summarise()`. The first warning was: ℹ In argument:
## `Age = chisq.test(Age, Orig$resp)$p.value`. Caused by warning in
## `chisq.test()`: ! Chi-squared approximation may be incorrect ℹ Run
## `dplyr::last_dplyr_warnings()` to see the 5 remaining warnings.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## # A tibble: 1 × 14
## resp Age Sex ChestPain RestBP Chol Fbs RestECG MaxHR ExAng
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.11e-65 0.0794 2.95e-6 1.18e-16 0.534 0.0951 1 0.00833 0.0718 9.51e-13
## # ℹ 4 more variables: Oldpeak <dbl>, Slope <dbl>, Ca <dbl>, Thal <dbl>
barplot(t(t(as.matrix(ChiSq.pval))), las=2, cex.names=1,
main="p-values from Chi-sq test of association")
abline(h=.05, col='red')
which(ChiSq.pval < .05) # Col num of variables w <.05 p-value
## [1] 1 3 4 8 10 11 12 13 14
list_cols[which(ChiSq.pval < .05)] # Top list of 'important' variables
## [1] "resp" "Sex" "ChestPain" "RestECG" "ExAng" "Oldpeak"
## [7] "Slope" "Ca" "Thal"
which(ChiSq.pval > .6) # Col num of variables w >.6 p-value
## [1] 7
list_cols[which(ChiSq.pval > .6)] # Top list of 'unimportatnt' variables
## [1] "Fbs"
CODE ONLY
(top)
##----------------------------------------
## 1. Load data from package MASS
library(tidyverse)
Heart <- read_csv("https://nmimoto.github.io/datasets/heart.csv")
Heart
class(Heart) # it's already tibble
print(Heart, n=100) # if you want to see more rows
print(Heart, width=1000) # if you want to see more cols
names(Heart)
# #0 "X1" index
# #1 "Age" age in years
# #2 "Sex" (Qualitative) sex (1 = male; 0 = female)
# #3 "ChestPain" (Qualitative) chest pain type
# -- Value 1: typical angina
# -- Value 2: atypical angina
# -- Value 3: non-anginal pain
# -- Value 4: asymptomatic
# #4 "RestBP" resting blood pressure (in mm Hg on admission to the hospital)
# #5 "Chol" serum cholestoral in mg/dl
# #6 "Fbs" (Qualitative) fasting blood sugar > 120 mg/dl (1 = true; 0 = false)
# #7 "RestECG" (Qualitative) resting electrocardiographic results
# -- Value 0: normal
# -- Value 1: having ST-T wave abnormality (T wave inversions and/or ST elevation or depression of > 0.05 mV)
# -- Value 2: showing probable or definite left ventricular hypertrophy by Estes' criteria
# #8 "MaxHR" maximum heart rate achieved
# #9 "ExAng" (Qualitative) exercise induced angina (1 = yes; 0 = no)
# #10 "Oldpeak" oldpeak = ST depression induced by exercise relative to rest
# #11 "Slope" the slope of the peak exercise ST segment
# #12 "Ca" number of major vessels (0-3) colored by flourosopy
# #13 "Thal" (Qualitative) Thalium stress test
# 3 = normal; 6 = fixed defect; 7 = reversable defect
# #14 "AHD" (Qualitative)(Response) Yes/No based on an angiographic test
# diagnosis of heart disease (angiographic disease status)
# -- Value 0: < 50% diameter narrowing
# -- Value 1: > 50% diameter narrowing
# (in any major vessel: attributes 59 through 68 are vessels)
##----------------------------------------
## 2. Data Preparation
### column "AHD" is the respoinse variable.
# table of the respoonse variable
table(Heart$AHD)
print(Heart, width=1000) # show all columns
# Rename medv column as resp.
Heart2 <- Heart %>%
select(-"index") %>% # remove col named "index"
rename(resp=AHD) %>% # rename the column
relocate(resp) # move "resp" to 1st column
print(Heart2, width=1000)
# Turn class of some columns into <factor>
Heart2 <- Heart2 %>%
mutate(resp = as.factor(resp), # Turn these columns to <factor> instead of <double>
Sex = as.factor(Sex),
ChestPain= as.factor(ChestPain),
Fbs = as.factor(Fbs),
RestECG = as.factor(RestECG),
ExAng = as.factor(ExAng),
Thal = as.factor(Thal))
print(Heart2, width=1000)
Heart2 %>% select(resp) %>% table
# note that level is "No/Yes", not "Yes/No". <---- Important!
##----------------------------------------
## 3. Routine Exploratory Analysis (class of resp should be <fct>)
Orig <- Heart2
resp.col.name <- "resp"
#- Check for N/A in data. Remove if there's any.
summary(Orig)
dim(Orig)
sum(is.na(Orig))
# If there is na in the data, run below
dim(Orig)
Orig <- Orig %>% na.omit()
dim(Orig)
##----------
##- Correlation Check
# Pick columns that are numeric
library(corrplot) # install.packages("corrplot")
Orig_num <- Orig %>% select_if(is.numeric)
cor(Orig_num)
corrplot::corrplot(cor(Orig_num))
corrplot::corrplot(cor(Orig_num), method="number")
##----------
##- Visualization
# pairs(Orig) equivalent
library(GGally) # install.packages("GGally")
GGally::ggpairs(Orig[, c(2:7, 1)], aes(color=resp, alpha=1))
GGally::ggpairs(Orig[, c(8:14, 1)], aes(color=resp, alpha=1))
###----------
###- Chi-sq test of association for each column
# Pick columns that has more than 1 unique value
list_cols <- Orig %>% summarise_all(function(x) length(unique(x))) %>%
gather() %>% filter(value>1) %>% pull(key)
# Apply chisq.test to those columns
ChiSq.pval <- Orig %>% select(list_cols) %>%
summarise_all(funs(chisq.test(., Orig$resp)$p.value))
ChiSq.pval <- Orig %>% summarise_all(funs(chisq.test(., Orig$resp)$p.value))
ChiSq.pval
barplot(t(t(as.matrix(ChiSq.pval))), las=2, cex.names=1,
main="p-values from Chi-sq test of association")
abline(h=.05, col='red')
which(ChiSq.pval < .05) # Col num of variables w <.05 p-value
list_cols[which(ChiSq.pval < .05)] # Top list of 'important' variables
which(ChiSq.pval > .6) # Col num of variables w >.6 p-value
list_cols[which(ChiSq.pval > .6)] # Top list of 'unimportatnt' variables