cardio-sem-echo-01.Rmd

---
title: "Análise Dataset de Cárdio"
author: "José Elvano Moraes"
date: "4/27/2021"
output:
  html_document: default
  pdf_document: default
---

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

```{r echo=FALSE}
suppressPackageStartupMessages(expr = library(tidyverse))
suppressPackageStartupMessages(expr = library(haven))
suppressPackageStartupMessages(expr = library(bnlearn))
suppressPackageStartupMessages(expr = library(Rgraphviz))
suppressPackageStartupMessages(expr = library(gRain))
suppressPackageStartupMessages(expr = library(Rmpfr))
suppressPackageStartupMessages(expr = library(readr))
suppressPackageStartupMessages(expr = library(precrec))
suppressPackageStartupMessages(expr = library(ROCR))
suppressPackageStartupMessages(expr = library(epiR))
```


```{r echo=FALSE}
suppressWarnings( 
  cardio <- read_csv("cardio_train.csv", 
                     col_types = cols(id = col_integer(), 
                     age = col_number(), 
                     gender = col_factor(levels = c("1", "2")), 
                     height = col_number(), 
                     weight = col_number(), 
                     ap_hi = col_number(), 
                     ap_lo = col_number(), 
                     cholesterol = col_factor(levels = c("1", "2", "3")), 
                     gluc = col_factor(levels = c("1", "2", "3")), 
                     smoke = col_factor(levels = c("0","1")), 
                     alco = col_factor(levels = c("0", "1")), 
                     active = col_factor(levels = c("0","1")), 
                     cardio = col_factor(levels = c("0","1")))
                     )
)
```

```{r echo=FALSE}
colunas <- c("age",         
             "gender",      
             "height",      
             "weight",      
             "ap_hi",       
             "ap_lo",       
             "cholesterol", 
             "gluc",        
             "smoke",      
             "alco",        
             "active",      
             "cardio") 
```

```{r echo=FALSE}
dd = filter(cardio, ap_hi>0, ap_hi>ap_lo, ap_lo > 0, ap_lo<300, ap_hi>0, ap_hi<300)
set.seed(2)


ccc <- select(dd, age, height, weight, ap_lo, ap_hi)
ccc2 <- discretize(ccc, method = 'hartemink', breaks = 4)
ddf <- cbind(select(dd, -age, -height, -weight, -ap_lo, -ap_hi), ccc2)
ddf <- na.omit(ddf)
ddf$id <- NULL

amostra = ddf[1:7000,]
amostra.test = ddf[7001:nrow(dd),]
amostra.test<-na.omit(amostra.test)
summary(ddf)
```

```{r echo=FALSE}
bn1 <- mmhc(ddf)
```

## Grafo mostrando como nodos coloridos as variáveis das quais a variavel **CARDIO** tem dependência estatística (*Markov blanket*) 
```{r echo=FALSE}
graphviz.plot(bn1, shape = 'rectangle', 
              main = "Markov Blanket - CARDIO",
              highlight = list(nodes=mb(bn1, 'cardio'), fill="lightblue"))
```

```{r echo=FALSE}
print("Fase de bootstrap...")
R = 100
arc.strength = boot.strength(ddf, R = R, algorithm = "mmhc")
rede.media = averaged.network(arc.strength)
```

## CDF da correlação (correlação de Pearson) entre os nodos do grafo
```{r echo=FALSE}
plot(arc.strength)
abline(v = 0.75, col = "tomato", lty = 2, lwd = 2)
abline(v = 0.85, col = "steelblue", lty = 2, lwd = 2)
```

## Grafo médio após 100 iterações 
```{r echo=FALSE}
par(mfrow = c(1, 2))
strength.plot(rede.media, 
              arc.strength, 
              shape = "rectangle")
```

```{r echo=FALSE}
# cado a rede nao seja totalmente direcional (nao seja um DAG) forcar-la a ser
rede.media.dagged <- pdag2dag(rede.media, 
                             ordering = names(ddf))
```

## Grafo médio tornado um *DAG*
```{r echo=FALSE}
graphviz.plot(rede.media.dagged, shape = 'rectangle', 
              main = "Markov Blanket - CARDIO",
              highlight = list(nodes=mb(rede.media.dagged, 'cardio'), fill="lightblue"))
```

## Distribuiçaõ estatística (PDF) de cada variável
```{r echo=FALSE}
fitted.m <- bn.fit(rede.media.dagged, amostra)

graphviz.chart(fitted.m,
               type = "barprob", 
               col = "darkblue",
               bg = "azure", 
               bar.col = "darkblue", 
               main = "DAG sem WL")
```


# predição
```{r echo=FALSE}
# rede.media.test.dagged <- pdag2dag(rede.media.test, 
#                                    ordering = names(amostra.test))

#USAR fit da rede da amostra learn, nao test
fitt <- bn.fit(rede.media.dagged, ddf)

amostra.test.resul <- amostra.test
amostra.test.resul$CARDIO.PRED <- -1
nodos <- amostra.test[, mb(rede.media.dagged, 'cardio')]
nodes.evidencias <-  mb(rede.media.dagged, 'cardio')

junction = compile(as.grain(fitt))
x<-c()
for (ii in c(1:nrow(nodos))) {
  jedu = setEvidence(propagate = TRUE, junction, 
                     nodes = nodes.evidencias, 
                     states = c(as.character(nodos[ii,1]), 
                                as.character(nodos[ii,2]), 
                                as.character(nodos[ii,3])#,
                                #as.character(nodos[ii,4]), 
                                #as.character(nodos[ii,5])
                                ))
  
  z <- querygrain(jedu, nodes = c("cardio"))
  zz <- z$cardio
  x[ii] <- zz[1]
  amostra.test.resul$CARDIO.PRED[ii] <- x[ii]
  
}
```

## Desempenho estatístico da rede bayesiana
```{r echo=FALSE}
amostra.test.resul$PREDICAO <- ifelse(amostra.test.resul$CARDIO.PRED>0.5, 0, 1)
v <- amostra.test.resul[,c('cardio', 'PREDICAO')]

# da lib epiR
print(epi.tests(table(v), conf.level = 0.95))
```
## Curva ROC
```{r echo=FALSE}
# Da library ROCR
pred <- prediction(1-x, v$cardio)
perf <- performance(pred,"tpr","fpr")
plot(perf,
     colorize=FALSE, 
     main="Curva ROC", 
     xlab = "Taxa de falsos positivos",
     ylab = "Taxa de verdeiro positivos")
```

## Outra implementação de curva  ROC e outras estatísticas
```{r echo=FALSE}
# da lib precrec
precrec_obj <- evalmod(scores = 1-x, labels = v$cardio)
autoplot(precrec_obj)

precrec_obj2 <- evalmod(scores = 1-x, labels = v$cardio, mode="basic")
autoplot(precrec_obj2)  
```

Þ