# Read datasets Cleveland_hd.csv into hd_data hd_data <- read.csv("datasets/Cleveland_hd.csv") # take a look at the first 5 rows of hd_data head(hd_data,5)

# load the tidyverse package library(dplyr) library(ggplot2) library(tidyr) # Use the 'mutate' function from dplyr to recode our data hd_data %>% mutate(hd = ifelse(class > 0, 1, 0))-> hd_data # recode sex using mutate function and save as hd_data hd_data %>% mutate(sex = factor(sex, levels = 0:1, labels = c('Female','Male')))-> hd_data;

# Does sex have an effect? Sex is a binary variable in this dataset, # so the appropriate test is chi-squared test hd_sex <- chisq.test(hd_data$sex,hd_data$hd) # Does age have an effect? Age is continuous, so we use a t-test hd_age <- t.test(hd_data$age~hd_data$hd) # What about thalach? Thalach is continuous, so we use a t-test hd_heartrate <- t.test(hd_data$thalach~hd_data$hd) # Print the results to see if p<0.05. print(hd_sex) print(hd_age) print(hd_heartrate)

# Recode hd to be labelled hd_data%>%mutate(hd_labelled = ifelse(hd == 0, 'No Disease', 'Disease')) -> hd_data # age vs hd ggplot(data = hd_data, aes(x = hd_labelled,y = age)) + geom_boxplot()

# sex vs hd ggplot(data = hd_data,aes(x=hd_labelled,fill=sex)) + geom_bar(position='fill')+ ylab("Sex %")

# max heart rate vs hd ggplot(data = hd_data,aes(x=hd_labelled,y=thalach)) + geom_boxplot()

# use glm function from base R and specify the family argument as binomial model <- glm(data = hd_data, hd~age+sex+thalach, family='binomial') # extract the model summary summary(model)

# load the broom package library('broom') # tidy up the coefficient table tidy_m <- tidy(model) tidy_m # calculate OR tidy_m$OR <- exp(tidy_m$estimate) # calculate 95% CI and save as lower CI and upper CI tidy_m$lower_CI <- exp(tidy_m$estimate - 1.96 * tidy_m$std.error) tidy_m$upper_CI <- exp(tidy_m$estimate + 1.96 * tidy_m$std.error) # display the updated coefficient table tidy_m

# get the predicted probability in our dataset using the predict() function pred_prob <- predict(model,hd_data, type = "response") # create a decision rule using probability 0.5 as cutoff and save the predicted decision into the main data frame hd_data$pred_hd <- ifelse(pred_prob>=0.5,1,0) # create a newdata data frame to save a new case information newdata <- data.frame(age = 45, sex = "Female", thalach = 150) # predict probability for this new case and print out the predicted value p_new <- predict(model,newdata, type = "response") p_new

# load Metrics package library(Metrics) # calculate auc, accuracy, clasification error auc <- auc(hd_data$hd,hd_data$pred_hd) accuracy <- accuracy(hd_data$hd,hd_data$pred_hd) classification_error <- ce(hd_data$hd,hd_data$pred_hd) # print out the metrics on to screen print(paste("AUC=", auc)) print(paste("Accuracy=", accuracy)) print(paste("Classification Error=", classification_error)) # confusion matrix table(hd_data$hd,hd_data$pred_hd, dnn=c('True Status','Predicted Status')) # confusion matrix