---
title: "Heart Disease - Data Analysis"
author: "Prof. Radhendushka Srivastava, Department of Mathematics, IIT Bombay"
date: "05/03/2022"
output:
  html_document: default
---

# Loading libraries
```{r,warning=F}
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(readxl))
suppressPackageStartupMessages(library(lessR))
```

# Reading data
```{r}
# Source: https://archive.ics.uci.edu/ml/datasets/heart+disease
heart <- read_xlsx("heart.xlsx")
heart$age <- as.numeric(heart$age)
heart$trestbps <- as.numeric(heart$trestbps)
heart$chol <- as.numeric(heart$chol)
heart$thalach <- as.numeric(heart$thalach)
heart$oldpeak <- as.numeric(heart$oldpeak)
heart$ca <- as.numeric(heart$ca)
```

# Discrete Variables
## 1) Univariate probablity mass function and visualization using bar plot & pie chart

### 1.1) P.M.F. of the variable "sex" (gender) and it's visualization
```{r}

# Spoken Tutorials:
# 15. Plotting Histograms and Pie Chart 
# 16. Plotting Bar Charts and Scatter Plot

prop.table(table(heart$sex))
# Bar plot
barplot(height = prop.table(table(heart$sex)), col = c("green","blue"), main = "Gender", xlab = "", ylab = "Proportion")
# Pie chart
pie(x = prop.table(table(heart$sex)), labels = names(prop.table(table(heart$sex))), col = c("green","blue"), radius = 1, main = "Gender")
```

### 1.2) P.M.F. of the variable "cp" (chest pain) and it's visualization
```{r}
prop.table(table(heart$cp))
# Bar plot
barplot(height = prop.table(table(heart$cp)), col = c("red","indianred1","violetred2","orangered3"), main = "Chest pain type", xlab = "", ylab = "Proportion")
# Pie chart
pie(x = prop.table(table(heart$cp)), labels = names(prop.table(table(heart$cp))), col = c("red","indianred1","violetred2","orangered3"), radius = 1, main = "Chest pain type")
```

### 1.3) P.M.F. of the variable "restecg" (resting ECG) and it's visualization
```{r}
prop.table(table(heart$restecg))
# Bar plot
barplot(height = prop.table(table(heart$restecg)), col = c("blue","slateblue1","steelblue2"), main = "Resting electrocardiographic results", xlab = "", ylab = "Proportion")
# Pie chart
pie(x = prop.table(table(heart$restecg)), labels = names(prop.table(table(heart$restecg))), col = c("blue","slateblue1","steelblue2"), radius = 1, main = "Resting electrocardiographic results")
```

### 1.4) P.M.F. of the variable "slope" (slope of the peak exercise ST segment) and it's visualization
```{r}
prop.table(table(heart$slope))
# Bar plot
barplot(height = prop.table(table(heart$slope)), col = c("green","limegreen","lawngreen"), main = "Slope of the peak exercise ST segment", xlab = "", ylab = "Proportion")
# Pie chart
pie(x = prop.table(table(heart$slope)), labels = names(prop.table(table(heart$slope))), col = c("green","limegreen","lawngreen"), radius = 1, main = "Slope of the peak exercise ST segment")
```

### 1.5) P.M.F. of the variable "thal" (thalasemia blood disorder) and it's visualization
```{r}
prop.table(table(heart$thal))
# Bar plot
barplot(height = prop.table(table(heart$thal)), col = c("peru","tan1","sandybrown"), main = "Thal", xlab = "", ylab = "Proportion")
# Pie chart
pie(x = prop.table(table(heart$thal)), labels = names(prop.table(table(heart$thal))), col = c("peru","tan1","sandybrown"), radius = 1, main = "Thal")
```

### 1.6) P.M.F. of variable "ca" (number of major vessels) and it's visualization
```{r}
prop.table(table(heart$ca))
# Bar plot
barplot(height = prop.table(table(heart$ca)), col = c("grey","lightcyan","ivory","honeydew"), main = "Proportion of major vessels", xlab = "", ylab = "Proportion")
# Pie chart
pie(x = prop.table(table(heart$ca)), labels = names(prop.table(table(heart$ca))), col = c("grey","lightcyan","ivory","honeydew"), radius = 1, main = "Proportion of major vessels")
```

## 2) Bivariate joint probability mass function and visualization using grouped bar plot
```{r}

# Spoken Tutorials:
# 23. Functions in R

# Function for cross tabulation (proportion)
cross_tab_proportion <- function(x,y)
{
  temp_table <- prop.table(table(x,y))
  temp_col <- colSums(temp_table)
  temp_row <- rowSums(temp_table)
  temp_sum <- sum(temp_table)
  temp <- as.data.frame(rbind(cbind(temp_table,temp_row),c(temp_col,temp_sum)))
  rownames(temp)[nrow(temp)] <- "total"
  colnames(temp)[ncol(temp)] <- "total"
  return(temp)
}
```

### 2.1) Bivariate joint p.m.f of gender and target along with grouped bar plot
```{r}
cross_tab_proportion(heart$sex,heart$target)
# Grouped bar plot
data <- prop.table(table(heart$target,heart$sex))
barplot(data, main = "Gender", xlab = "", ylab = "Proportion",
        col = c("blue", "red", "blue", "red"), beside = TRUE)
legend("topleft", rownames(data), fill = c("blue", "red"), xpd=TRUE, inset=c(0, -.1), horiz = T, bty = "n")
```

### 2.2) Bivariate joint p.m.f of chest pain and target along with grouped bar plot
```{r}
cross_tab_proportion(heart$cp,heart$target)
# Grouped bar plot
data <- prop.table(table(heart$target,heart$cp))
barplot(data, main = "Chest pain type", xlab = "", ylab = "Proportion",
        col = c("blue", "red", "blue", "red"), beside = TRUE)
legend("topleft", rownames(data), fill = c("blue", "red"), xpd=TRUE, inset=c(0, -.1), horiz = T, bty = "n")
```

### 2.3) Bivariate joint p.m.f of resting ECG and target along with grouped bar plot
```{r}
cross_tab_proportion(heart$restecg,heart$target)
# Grouped bar plot
data <- prop.table(table(heart$target,heart$restecg))
barplot(data, main = "Resting electrocardiographic results", xlab = "", ylab = "Proportion",
        col = c("blue", "red", "blue", "red"), beside = TRUE)
legend("topleft", rownames(data), fill = c("blue", "red"), xpd=TRUE, inset=c(0, -.1), horiz = T, bty = "n")
```

### 2.4) Bivariate joint p.m.f of the slope of the peak exercise ST segment and target along with grouped bar plot
```{r}
cross_tab_proportion(heart$slope,heart$target)
# Grouped bar plot
data <- prop.table(table(heart$target,heart$slope))
barplot(data, main = "Slope of the peak exercise ST segment", xlab = "", ylab = "Proportion",
        col = c("blue", "red", "blue", "red"), beside = TRUE)
legend("topleft", rownames(data), fill = c("blue", "red"), xpd=TRUE, inset=c(0, -.1), horiz = T, bty = "n")
```

### 2.5) Bivariate joint p.m.f of thalesemia blood disorder and target along with grouped bar plot
```{r}
cross_tab_proportion(heart$thal,heart$target)
# Grouped bar plot
data <- prop.table(table(heart$target,heart$thal))
barplot(data, main = "Thal", xlab = "", ylab = "Proportion",
        col = c("blue", "red", "blue", "red"), beside = TRUE)
legend("topleft", rownames(data), fill = c("blue", "red"), xpd=TRUE, inset=c(0, -.1), horiz = T, bty = "n")
```

### 2.6) Bivariate joint p.m.f of the number of major vessels and target along with grouped bar plot
```{r}
cross_tab_proportion(heart$ca,heart$target)
# Grouped bar plot
data <- prop.table(table(heart$target,heart$ca))
barplot(data, main = "Proportion of major vessels", xlab = "", ylab = "Proportion",
        col = c("blue", "red", "blue", "red"), beside = TRUE)
legend("topright", rownames(data), fill = c("blue", "red"), xpd=TRUE, inset=c(0, -.1), horiz = T, bty = "n")
```

## 3) $\chi^2$ test: Testing association between "target" and other discrete variables 
```{r}
# Function for cross tabulation (count)
cross_tab_count <- function(x,y)
{
  temp_table <- table(x,y)
  temp_col <- colSums(temp_table)
  temp_row <- rowSums(temp_table)
  temp_sum <- sum(temp_table)
  temp <- as.data.frame(rbind(cbind(temp_table,temp_row),c(temp_col,temp_sum)))
  rownames(temp)[nrow(temp)] <- "total"
  colnames(temp)[ncol(temp)] <- "total"
  return(temp)
}
```

### 3.1) Testing association between gender and target
```{r}
cross_tab_count(heart$sex,heart$target)
# Chi-squared test
chisq.test(heart$sex,heart$target,correct = F)
```

### 3.2) Testing association between chest pain and target
```{r}
cross_tab_count(heart$cp,heart$target)
# Chi-squared test
chisq.test(heart$cp,heart$target)
```

### 3.3) Testing association between resting ECG and target
```{r}
cross_tab_count(heart$restecg,heart$target)
# Chi-squared tests
chisq.test(heart$restecg,heart$target)
```

### 3.4) Testing association between the slope of the peak exercise ST segment and target
```{r}
cross_tab_count(heart$slope,heart$target)
# Chi-squared test
chisq.test(heart$target,heart$slope)
```

### 3.5) Testing association between thalsemia blood disorder and target
```{r}
cross_tab_count(heart$thal,heart$target)
# Chi-squared test
chisq.test(heart$thal,heart$target)
```

### 3.6) Testing association between the number of major vessels and target
```{r}
cross_tab_count(heart$ca,heart$target)
# Chi-squared test
chisq.test(heart$ca,heart$target)
```

## 4) Testing of equality of two proportions
### 4.1) Testing whether the proportion of sick male and sick female are equal or not
```{r}

# Spoken Tutorials:
# 19. Data Manipulation using dplyr Package
# 20. More Functions in dplyr Package
# 21. Pipe Operator

sick_female <- as.vector(table(heart %>% filter(target=="sick") %>% select(sex)))[1]
sick_male <- as.vector(table(heart %>% filter(target=="sick") %>% select(sex)))[2]
total_female <- length(unlist(heart %>% select(sex) %>% filter(sex=="fem")))
total_male <- length(unlist(heart %>% select(sex) %>% filter(sex=="male")))
proportion_test <- function(n1,N1,n2,N2)
{
  p1 <- n1/N1
  p2 <- n2/N2
  p  <- (n1+n2)/(N1+N2)
  se <- sqrt(p*(1-p)*((1/N1)+1/(N2)))
  result <- (p1-p2)/se
  return(cat("p1:",p1,"p2:",p2,"z-statistic:",result,"p-value:",c(2*(1-pnorm(abs(result))))))
}
proportion_test(sick_female,total_female,sick_male,total_male)
```

### 4.2) Testing equality of two proportions by using a built-in R function
```{r}
sick_female <- as.vector(table(heart %>% filter(target=="sick") %>% select(sex)))[1]
sick_male <- as.vector(table(heart %>% filter(target=="sick") %>% select(sex)))[2]
total_female <- length(unlist(heart %>% select(sex) %>% filter(sex=="fem")))
total_male <- length(unlist(heart %>% select(sex) %>% filter(sex=="male")))
prop.test(x=c(sick_female,sick_male), n=c(total_female,total_male),alternative = "two.sided", correct = FALSE)
```

# Continuous variables
## 5) Basic statistics of continuous variables
```{r}
# Mean, standard deviation, minimum, maximum, range, quartiles, IQR, skewness and kurtosis
statistics <- function(x)
{
  mean <- mean(x)
  standard_deviation <- sd(x)
  min <- min(x)
  max <- max(x)
  range <- max - min
  first_quartile <- quantile(x)[2]
  third_quartile <- quantile(x)[4]
  IQR <- third_quartile - first_quartile
  skewness <- mean((x-mean(x))^3)/(sd(x)^3)
  kurtosis <- (mean((x-mean(x))^4))/sd(x)^4
  
  return(cat("\t","Mean: ",mean,"\n\t","Standard deviation: ",standard_deviation,"\n\t","Minimum value: ",min,"\n\t","Maximum value: ",max,"\n\t","Range: ",range,"\n\t","First Quartile: ",first_quartile,"\n\t","Third Quartile: ",third_quartile,"\n\t","Interquartile Range: ",IQR,"\n\t","Skewness: ",skewness,"\n\t","Kurtosis: ",kurtosis))
}
```

### 5.1) Basic statistics of the variable "age"
```{r}
statistics(heart$age)
```

### 5.2) Basic statistics of the variable "chol" (Cholestrol)
```{r}
statistics(heart$chol)
```

### 5.3) Basic statistics of the variable "thalach" (Maximum heart rate achieved)
```{r}
statistics(heart$thalach)
```

### 5.4) Basic statistics of the variable "oldpeak" (ST depression induced by exercise relative to rest)
```{r} 
statistics(heart$oldpeak)
```

## 6) Testing equality of mean of two population
### 6.1) Testing equality of the mean of "chol" among sick and healty group
```{r}
# Test of equality of mean
chol_sick <- unlist(heart %>% filter(target=="sick") %>% select(chol))
chol_healthy <- unlist(heart %>% filter(target=="healthy") %>% select(chol))
mean_test <- function(x1,x2)
{
  x_bar_1 <- mean(x1)
  x_bar_2 <- mean(x2)
  n1=length(x1)
  n2=length(x2)
  sd_1 <- sd(x1)
  sd_2 <- sd(x2)
  se <- sqrt(((n1-1)*c(sd_1^2)+(n2-1)*c(sd_2^2))/(n1+n2-2))*sqrt(c(1/n1)+c(1/n2))
  t_stat <- (x_bar_1-x_bar_2)/se
  return(cat("x_bar 1:",x_bar_1,"x_bar 2:",x_bar_2,"t-statistic:",t_stat,"se:", se,"df:",c(n1+n2-2),"p-value:",2*c(1-pt(abs(t_stat),df=c(n1+n2-2)))))
}
mean_test(chol_sick,chol_healthy)

```

### 6.2) Testing equality of mean using a built-in R function
```{r}
chol_sick <- unlist(heart %>% filter(target=="sick") %>% select(chol))
chol_healthy <- unlist(heart %>% filter(target=="healthy") %>% select(chol))
A=t.test(chol_sick,chol_healthy,var.equal=T)
A
```