heart_disease.Rmd

---
title: "project"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



# log normal distribution


```{r}
heart_disease<-read.csv('./heart.csv')
```
```{r}
heart_disease$HeartDisease<-as.factor(heart_disease$HeartDisease)
heart_disease$ChestPainType<-as.factor(heart_disease$ChestPainType)
heart_disease$ExerciseAngina<-as.factor(heart_disease$ExerciseAngina)
heart_disease$ST_Slope<-as.factor(heart_disease$ST_Slope)
heart_disease$Sex<-as.factor(heart_disease$Sex)
summary(heart_disease)
```


```{r}
head(heart_disease)
```
```{r}
library(ggplot2)
ggplot(data = heart_disease,  aes(x = Age,  y = RestingBP, color =  
                 Sex)) + 
     geom_point()
```
```{r}
library(GGally)
ggpairs(heart_disease)
```


```{r}
head(heart_disease)
```

```{r}
heart_logit<-glm(HeartDisease~., data = heart_disease, family = binomial)
summary(heart_logit)
```
```{r}
heart_logit_reduced = step(heart_logit)
```

```{r}
summary(heart_logit_reduced)
```

```{r}
LLR = 2 * (logLik(heart_logit) - logLik(heart_logit_reduced))
cat(logLik(heart_logit), logLik(heart_logit_reduced))
df = heart_logit_reduced$df.residual - heart_logit$df.residual
p_value = 1 - pchisq(LLR, df)
print(df)
cat('The loglikelihood ratio test statistic is', LLR, 'and the p-value is ', p_value, 'with df', df)
```

```{r}

heart_logit_all_square = glm(HeartDisease~.  +I(MaxHR^2) +I(Age^2) + I(Oldpeak^2)+ I(RestingBP^2) +I(Cholesterol^2), family = binomial, data = heart_disease)
summary(heart_logit_all_square)
```
```{r}
heart_logit_square_reduced = step(heart_logit_all_square)
summary(heart_logit_square_reduced)
```
```{r}
LLR = 2 * (logLik(heart_logit_all_square) - logLik(heart_logit_square_reduced))
df = heart_logit_square_reduced$df.residual - heart_logit_all_square$df.residual
p_value = 1 - pchisq(LLR, df)
print(logLik(heart_logit_all_square))
print( logLik(heart_logit_square_reduced))
cat('The loglikelihood ratio test statistic is', LLR, 'and the p-value is ', p_value)
```

```{r}
heart_logit_all_interaction = glm(HeartDisease~. +(MaxHR + Age + Oldpeak + RestingBP)^2 +I(MaxHR^2) +I(Age^2) + I(Oldpeak^2)+ I(RestingBP^2), family = binomial, data = heart_disease)
summary(heart_logit_all_interaction)
```
```{r}
heart_logit_interaction_reduced = step(heart_logit_all_interaction, direction = "both")
```
```{r}
summary(heart_logit_interaction_reduced)
```
```{r}
heart_logit_interaction_reduced = glm(formula = HeartDisease ~ Age + Sex + ChestPainType + Cholesterol + 
    FastingBS + MaxHR + ExerciseAngina + Oldpeak + ST_Slope + 
    I(Oldpeak^2) + MaxHR:Oldpeak, family = binomial, data = heart_disease)
summary(heart_logit_interaction_reduced)
```



```{r}
LLR = 2 * (logLik(heart_logit_all_interaction) - logLik(heart_logit_interaction_reduced))
df = heart_logit_interaction_reduced$df.residual - heart_logit_all_interaction$df.residual
p_value = 1 - pchisq(LLR, df)
print(logLik(heart_logit_all_interaction))
print(logLik(heart_logit_interaction_reduced))
print(df)
cat('The loglikelihood ratio test statistic is', LLR, 'and the p-value is ', p_value)
```

```{r}
null_model = glm(HeartDisease~ 1, data = heart_disease, family = binomial)

Pse_R2 = (deviance(null_model) - deviance(heart_logit_square_reduced)) /deviance(null_model)
cat('Pseudo-R square:', Pse_R2, '\n')
```
```{r}
null_model = glm(HeartDisease~ 1, data = heart_disease, family = binomial)

Pse_R2 = (deviance(null_model) - deviance(heart_logit_interaction_reduced)) /deviance(null_model)
cat('Pseudo-R square:', Pse_R2, '\n')
```
```{r}

null_model = glm(HeartDisease~ 1, data = heart_disease, family = binomial)

Pse_R2 = (deviance(null_model) - deviance(heart_logit_reduced)) /deviance(null_model)
cat('Pseudo-R square:', Pse_R2, '\n')
```
```{r}
library(pROC)
invisible(plot(roc(heart_disease$HeartDisease,
                   fitted(heart_logit_reduced)),
               col = "#5da492", 
               main = "ROC curves:3 models",
               legend = 'Logistic model of pure covariates'))

invisible(plot(roc(heart_disease$HeartDisease,
                   fitted(heart_logit_square_reduced)),
               print.auc = T, 
               col = "#206376", 
               add = T))
invisible(plot(roc(heart_disease$HeartDisease,
                   fitted(heart_logit_interaction_reduced)),
               col = "#bee7a3", 
               add = T))
color = c("#5da492","#206376","#bee7a3")
```

```{r}
roc1 = roc(heart_disease$HeartDisease,
                   fitted(heart_logit_reduced))
roc2 = roc(heart_disease$HeartDisease,
                   fitted(heart_logit_square_reduced))
roc3 = roc(heart_disease$HeartDisease,
                   fitted(heart_logit_interaction_reduced))



```
```{r}
library(reshape2)
library(plotROC)
ROC= data.frame(h = heart_disease$HeartDisease, model1 =fitted(heart_logit_reduced), model2=fitted(heart_logit_square_reduced),model3=fitted(heart_logit_interaction_reduced) )
ROC$h = as.numeric(as.character(ROC$h))
longtest <- melt_roc(ROC, "h", c("model3", "model2", "model1"))
head(longtest)
ggplot(longtest, aes(d = D, m = M, color = name), main = "3 models ROC curve", xlab = "False Positive Rate (1-Specificity)", ylab = "True Positive Rate (Sensitivity)") + 
    geom_roc() +
    style_roc() +  
    ggtitle('ROC Plot for Three Models')
   
```


### Deviance test
```{r}
s = summary(heart_logit_square_reduced)
dev = deviance(heart_logit_square_reduced)

p_value = 1-pchisq(dev, s$df.residual)
cat('The deviance of the model is : ', dev,'with degree freedom of ', s$df.residual, '\n')
cat('p_value is :', p_value)
```

```{r}
s = summary(heart_logit_interaction_reduced)
dev = deviance(s)

p_value = 1-pchisq(dev, s$df.residual)
cat('The deviance of the model is : ', dev,'with degree freedom of ', s$df.residual, '\n')
cat('p_value is :', p_value)
```


```{r}
s = summary(heart_logit_reduced)
dev = deviance(s)

p_value = 1-pchisq(dev, s$df.residual)
cat('The deviance of the model is : ', dev,'with degree freedom ofs ', s$df.residual, '\n')
cat('p_value is :', p_value)
```

```{r}    
library(arm)
b1 = binnedplot(fitted(heart_logit_reduced), 
        residuals(heart_logit_reduced, type = "response"), 
        nclass = NULL, 
        xlab = "Expected Values", 
        ylab = "Average residual", 
        main = "Binned residual plot for model 1", 
        cex.pts = 0.8, 
        col.pts = "red", 
        col.int = "blue",
    )
b2 = binnedplot(fitted(heart_logit_square_reduced), 
        residuals(heart_logit_square_reduced, type = "response"), 
        nclass = NULL, 
        xlab = "Expected Values", 
        ylab = "Average residual", 
        main = "Binned residual plot for model 2", 
        cex.pts = 0.8, 
        col.pts = "red", 
        col.int = "blue",
    )
b3 = binnedplot(fitted(heart_logit_interaction_reduced), 
        residuals(heart_logit_interaction_reduced, type = "response"), 
        nclass = NULL, 
        xlab = "Expected Values", 
        ylab = "Average residual", 
        main = "Binned residual plot for model 3", 
        cex.pts = 0.8, 
        col.pts = "red", 
        col.int = "blue",
    )
color = c("#5da492","#206376","#bee7a3")
        
```
```{r}
library(ggplot2)
ggpubr::ggarrange(b1, b2, b3, 
                  labels = c("Model 1", "Model 2", "Model 3"),
                  ncol = 1, nrow = 3)
```

```{r}

plot(heart_logit_square_reduced)
```



```{r}
heart_logit_square_reduced$coefficients
confint.default(heart_logit_square_reduced)
```


```{r}
exp(heart_logit_square_reduced$coefficients)
exp(confint.default(heart_logit_square_reduced))
```