rm(list=ls())
##
## Teorema de Bayes 
##   - probabilidades calculadas em dados
##
## Cálculos dos exemplos do teste de diagnóstico
##

## Valores necessários (neste caso fixados)
prev <- 0.02  ## prevalência
sens <- 0.90  ## sensibilidade 
esp  <- 0.80  ## especificidade

##
## Cálculos dos Resultados dos slides
## 
pr <- c(S = 1-prev, D = prev)
Tmat <- rbind(
    S = c(N = esp, P = 1-esp),
    D = c( N = 1-sens, P = sens)
)
pr          ## informação "a priori"
Tmat        ## características dos dados (verossimilhança)
pr*Tmat     ## combinando (conjunta)
## cenários para prevalência de 2%
post <- prop.table(pr*Tmat, mar = 2) 
post        ## informação pós teste (posteriori)
post[,2]    ## selecionando para resultado do teste positivo:
post[,1]    ## selecionando para resultado do teste negativo:
## ou selecionado por nome
post[, "P"]
post[, "N"]

## obtendo dois positivos em testes consecutivos 
(pos1 <- prop.table(pr*Tmat, mar = 2)[,"P"])
(pos2 <- prop.table(pos1*Tmat, mar = 2)[,"P"])

## e se forem tres positivos ...
(pos1 <- prop.table(pr*Tmat, mar = 2)[,"P"])
(pos2 <- prop.table(pos1*Tmat, mar = 2)[,"P"])
(pos3 <- prop.table(pos2*Tmat, mar = 2))
pos3[,"P"]
##
## Teorema de Bayes 
##   - probabilidades calculadas em dados
##
## Exemplo 1: classificando com uma preditora
##
df1 <- read.csv("http://www.leg.ufpr.br/~paulojus/CE315/dados/TesteDiag-dados1.csv")
str(df1)
dim(df1)
head(df1)

## prevalência: P[D]
with(df1, prop.table(table(Classe)))

## sensibilidade: P[+|D]
with(df1, prop.table(table(Classe, X1), mar=1))
with(subset(df1, Classe == "Doente"), prop.table(table(X1)))

## especificidade: P[-|S]
with(df1, prop.table(table(Classe, X1), mar=1))
with(subset(df1, Classe == "Sadio"), prop.table(table(X1)))

## vpp: P[D|+]
with(df1, prop.table(table(Classe, X1), mar=2))
with(subset(df1, X1 == "Positivo"), prop.table(table(Classe)))

## vpn: P[S|-]
with(df1, prop.table(table(Classe, X1), mar=2))
with(subset(df1, X1 == "Positivo"), prop.table(table(Classe)))

rm(list=ls())
##
## Exemplo 2: classificando com três preditoras
##
df <- read.csv("http://www.leg.ufpr.br/~paulojus/CE315/dados/TesteDiag-dados3.csv")
dim(df)
head(df)

## priori: prevalência
with(df, prop.table(table(Classe)))

## verossimilhança: sens., esp. e complementares como uma tabela - matrix
with(df, prop.table(table(Classe, X1), mar = 1))
with(df, prop.table(table(Classe, X2), mar = 1))
with(df, prop.table(table(Classe, X3), mar = 1))

## posterioris para cada dado (teste) separadamente
with(subset(df, X1 == "Positivo"), prop.table(table(Classe)))
with(subset(df, X1 == "Negativo"), prop.table(table(Classe)))

with(df, prop.table(table(Classe, X1), mar = 2))
with(df, prop.table(table(Classe, X2), mar = 2))
with(df, prop.table(table(Classe, X3), mar = 2))

## posteriori após os 3 testes
res <- c("Positivo", "Positivo", "Positivo")

## Forma 1:
## exata na base
df.res <- subset(df, X1 == res[1] & X2 == res[2] & X3 == res[3])
dim(df.res)
with(df.res, prop.table(table(Classe)))

## Forma 2:
## estimada supondo independência (condicional):

## calculando dados novamente
(pr <- with(df, as.vector(prop.table(table(Classe)))))
(tb1 <- with(df, prop.table(table(Classe, X1), mar = 1)))
(tb2 <- with(df, prop.table(table(Classe, X2), mar = 1)))
(tb3 <- with(df, prop.table(table(Classe, X3), mar = 1)))

## ideia para posteriori
with(df, prop.table(table(Classe, X1), mar = 2))
with(df, prop.table(pr*tb1, mar = 2))

## posterioris sequenciais (produto)
(post1 <- with(df, prop.table(pr*tb1, mar = 2))[,res[1]])
(post2 <- with(df, prop.table(post1*tb2, mar = 2))[,res[2]])
(post3 <- with(df, prop.table(post2*tb3, mar = 2))[,res[3]])

## Comparar com resultado de Forma 1

## Juntando tudo para calcular os todos possíveis perfis de resposta

post.fun <- function(res, df){
    ## res é um vetor com os resultados dos testes
    pr <- with(df, as.vector(prop.table(table(Classe))))
    tb1 <- with(df, prop.table(table(Classe, X1), mar = 1))
    tb2 <- with(df, prop.table(table(Classe, X2), mar = 1))
    tb3 <- with(df, prop.table(table(Classe, X3), mar = 1))
    
    post1 <- with(df, prop.table(pr*tb1, mar = 2))[,res[1]]
    post2 <- with(df, prop.table(post1*tb2, mar = 2))[,res[2]]
    post3 <- with(df, prop.table(post2*tb3, mar = 2))[,res[3]]
    
    post3
}

res.df <- expand.grid(
    X1 = c("Negativo", "Positivo"),
    X2 = c("Negativo", "Positivo"),
    X3 = c("Negativo", "Positivo")
)
res.df

apply(res.df, 1, post.fun, df = df)
cbind(res.df, round(t(apply(res.df, 1, post.fun, df = df)), dig=4))
##
## Exemplo 3: separando conjuntos de treino e teste
##
rm(list=ls())
df <- read.csv("http://www.leg.ufpr.br/~paulojus/CE315/dados/TesteDiag-dados3.csv")
dim(df)
head(df)

set.seed(2025)
IND <- sample(1:nrow(df), round(0.3*nrow(df))) ## para separar 30%
df.teste <- df[IND,]
df.treino <- df[-IND,]
dim(df.treino)
dim(df.teste)

## calculando dados novamente agora somente para o conjunto de treino
(pr <- with(df.treino, as.vector(prop.table(table(Classe)))))
(tb1 <- with(df.treino, prop.table(table(Classe, X1), mar = 1)))
(tb2 <- with(df.treino, prop.table(table(Classe, X2), mar = 1)))
(tb3 <- with(df.treino, prop.table(table(Classe, X3), mar = 1)))

## Predição com conjunto de teste
head(df.teste)
## posteriori após os 3 testes
## posterioris sequenciais (produto)
pred.f <- function(res){
    post1 <- prop.table(pr*tb1, mar = 2)[,res[1]]
    post2 <- prop.table(post1*tb2, mar = 2)[,res[2]]
    post3 <- prop.table(post2*tb3, mar = 2)[,res[3]]
    return(post3)
} 
pred.f(c("Positivo", "Positivo", "Positivo")) ## qual a classificação?
pred.f(c("Negativo", "Positivo", "Positivo")) ## qual a classificação?
pred.f(c("Positivo", "Negativo", "Positivo")) ## qual a classificação?
## por que os dois últimos não deram o mesmo resultado?
## quando dariam o mesmo resultado?

## agora todas do conjunto de teste
df.teste$probs <- t(apply(df.teste[,c("X1","X2","X3")], 1, pred.f))
head(df.teste$probs)
df.teste$predClasse <- ifelse(df.teste$probs[,1] > 0.5, "Doente", "Sadio")
head(df.teste, 20)
##
## Avaliando predição
with(df.teste, table(Classe, predClasse))
## pode-se mudar o "ponto de corte"
df.teste$predClasse <- ifelse(df.teste$probs[,1] > 0.25, "Doente", "Sadio")
with(df.teste, table(Classe, predClasse))
df.teste$predClasse <- ifelse(df.teste$probs[,1] > 0.5, "Doente", "Sadio")

MC <- with(df.teste, table(Classe, predClasse))
MC
sum(diag(MC))/sum(MC) ## acurácia
MC[2,2]/sum(MC[,2]) ## VPP
MC[1,1]/sum(MC[,1]) ## VPN
MC[2,2]/sum(MC[2,]) ## sensibilidade
MC[1,1]/sum(MC[1,]) ## especificidade
##
## Curva ROC
##
require(pROC)
roc1 <- with(df.teste,
             roc(Classe, probs[,1],
                 levels = c("Sadio", "Doente"),
                 direction = "<"))
roc1
plot(roc1, col="blue", lwd=2)
abline(a=0, b=1, lty=2, col="gray")
auc(roc1)
##

##
## Mas isto estas avaliações não fazem muito sentido
## neste exemplo com poucas combinações de preditores
##
table(df.teste$probs)  ## apenas 8x2=16 valores possíveis

res.df <- expand.grid(
    X1 = c("Negativo", "Positivo"),
    X2 = c("Negativo", "Positivo"),
    X3 = c("Negativo", "Positivo")
)
res.df
apply(res.df, 1, pred.f)
##
## Utilizando um pacote para Naïve Bayes
##
require(naivebayes)
nb <- naive_bayes(Classe ~ X1 + X2 + X3, data = df.treino)
nb
df.teste$nbprob <- predict(nb, newdata = df.teste, type = "prob")
df.teste$nbpred <- predict(nb, newdata = df.teste, type = "class")
## ou fazendo simplesmente para as 8 classes possíveis
predict(nb, newdata = res.df, type = "prob")
##
## Conferindo resultados do meu código com o do pacotes
##
head(df.teste, 20)[, c("probs", "nbprob")]
with(df.teste, table(predClasse, nbpred))
##
##