#################################################
#######   AULA 23/04                      #######
#######   MAE5776 - IME - 2o.Sem-2020     #######    
#######   Ana Gabriela P. de Vasconcelos  #######
#################################################

library(clusterGeneration)
library(biotools)
library(MASS)
library(dplyr)
library(MVN)

####### Análise discriminante #######

###### Análise discriminante de Fisher ######
#Suposições
#População estratificada em G grupos
#Independências
#Homocedasticidade
#Redução de dimensionalidade
#m < min(n,p,G-1)

#Exemplo questão 2 Lista 2
#-- Cenário 1: variâncias com sinais iguais --#

gera.dados <- function(G,ng,cenario){
  #Matriz SigmaB
  set.seed(541); varianciasB <- runif(2, min=3000, max=5000)
  set.seed(53);corrB <- rcorrmatrix(2, alphad=1)
  B <- diag(sqrt(varianciasB)) %*% corrB %*% t(diag(sqrt(varianciasB)))
  
  #Matriz SigmaW
  set.seed(313); varianciasW <- runif(2, min=20, max=200)
  set.seed(2190);corrW <- rcorrmatrix(2, alphad=1);corrW
  if(cenario == 1) corrW <- abs(corrW)
  W <- diag(sqrt(varianciasW)) %*% corrW %*% t(diag(sqrt(varianciasW)));W
  
  set.seed(542); mu <- mvrnorm(n=G, mu=c(0,0), Sigma = B)
  dados <- data.frame(Grupo=NULL, Y1=NULL,Y2=NULL)
  for(i in 1:G){
    set.seed(i); normal <- mvrnorm(n=ng, mu = mu[i,], Sigma = W)
    dados <- rbind(dados,cbind(rep(i,ng),normal))
  }
  colnames(dados) <- c("Grupo","Y1","Y2")
  dados$Grupo <- as.factor(dados$Grupo)
  return(dados)
}

#--- G=2 ---#
G1 <- 2
ng <- 100
dados1 <- gera.dados(G=G1,ng=ng,cenario=1)
plot(dados1$Y1, dados1$Y2, col=dados1$Grupo)

#Teste de homogeneidade de variâncias
boxM(dados1[,2:3],dados1$Grupo)  #Hipótese de homocedasticidade aceita

#Calculando (slide 13 aula 10)
medias <- aggregate(cbind(Y1,Y2) ~ Grupo, dados1, FUN=mean)
variancia1 <- cov(dados1[dados1$Grupo == '1',-1])
variancia2 <- cov(dados1[dados1$Grupo == '2',-1])
Sc <- ( (ng - 1) * variancia1 + (ng - 1) * variancia2 )/ (ng + ng -2)

#min(n,p,G-1) = 1
cargas <- as.matrix(solve(Sc)) %*% t(as.matrix(medias[1,-1] - medias[2,-1]))
escores <- as.matrix(dados1[,-1])  %*% cargas

regra.c <- c((1/2) * (as.matrix(medias[1,-1] - medias[2,-1])) %*% solve(Sc)  
             %*% t(as.matrix(medias[1,-1] + medias[2,-1]))); regra.c
class <- ifelse(escores[,1] > regra.c, 1,2)
table(dados1[,1],class)

plot(escores, col = class)

#Na prática não é necessário calcular na mão
fit <- lda(Grupo ~ Y1 + Y2, data = dados1, prior = c(1,1)/2)

fit$scaling
cargas #diferentes padronizações (slide 25 aula 10?)


fit.escores <- predict(fit, dados1[,-1])
names(fit.escores)
table(fit.escores$class,dados1$Grupo)

#--- G > 2 ---#
G2 <- 5
dados2 <- gera.dados(G=G2,ng=ng,cenario=1)
plot(dados2[,-1], col = dados2$Grupo)

fit.manova <- manova(cbind(Y1,Y2) ~ Grupo, data = dados2)

Sb <- summary(fit.manova)$SS$Grupo
Sw <- summary(fit.manova)$SS$Residual

fit.eigen <- eigen(solve(Sw)%*%Sb)
cargas2 <- fit.eigen$vectors
escores2 <- as.matrix(dados2[,-1])  %*% cargas2
plot(escores2, col = dados2[,1])

#Classificação (página 628 Johnson)
mu.dados2 <- aggregate(cbind(Y1,Y2)~Grupo, data = dados2, FUN = mean)
mu.escores <- as.data.frame(as.matrix(mu.dados2[,-1]) %*% t(cargas2))
mu.escores$Grupo <- c(1:G2)
colnames(mu.escores) <- c("Y1","Y2","Grupo")

class2 <- rep(NA, nrow(escores2))
for(i in c(1:nrow(escores2))){ #Função simples que não considera casos de distancias iguais
  distancias <- as.vector(dist(rbind(escores2[i,],mu.escores[,-3])))[c(1:G2)]
  class2[i] <- which(distancias == min(distancias))
}
plot(escores2, col = class2)

#min(n,p,G-1) = 2
fit2 <- lda(Grupo ~ Y1 + Y2, data = dados2, prior = rep(1,G2)/G2)
fit2$scaling

fit2.escores <- predict(fit2,dados2[,-1])
ct2 <- table(fit2.escores$class, dados2$Grupo);ct2
sum(diag(prop.table(ct2))) # prop. classificação correta total

plot(fit2.escores$x, col = fit2.escores$class)
plot(fit2.escores$x, col = dados2$Grupo)

#Comparação entre cenários
dados.cen1 <- gera.dados(G=G2,ng=ng,cenario=1)
dados.cen2 <- gera.dados(G=G2,ng=ng,cenario=2)

par(mfrow=c(1,2))
plot(dados.cen1[,-1], col = dados.cen1$Grupo, main = "Mesmo sinal")
plot(dados.cen2[,-1], col = dados.cen2$Grupo, main = "Sinais Opostos")
par(mfrow=c(1,1))

#Cargas
fit.cen1 <- lda(Grupo ~ Y1 + Y2, data = dados.cen1, prior = rep(1,G2)/G2)
fit.cen1$scaling 
fit.cen1$svd

fit.cen2 <- lda(Grupo ~ Y1 + Y2, data = dados.cen2, prior = rep(1,G2)/G2)
fit.cen2$scaling 
fit.cen2$svd


#Escores
fit.cen1.escores <- predict(fit.cen1, dados.cen1[,-1])
ct.cen1 <- table(fit.cen1.escores$class, dados.cen1$Grupo);ct.cen1
sum(diag(prop.table(ct.cen1))) # prop. classificação correta total

fit.cen2.escores <- predict(fit.cen2, dados.cen2[,-1])
ct.cen2 <- table(fit.cen2.escores$class, dados.cen2$Grupo);ct.cen2
sum(diag(prop.table(ct.cen2))) # prop. classificação correta total


###### Análise discriminante Solução Probabilística ######

#https://www.kaggle.com/imnikhilanand/heart-attack-prediction
#age: age in years
#trestbps: resting blood pressure (in mm Hg on admission to the hospital)
#chol: serum cholestoral in mg/dl
#thalach: maximum heart rate achieved
#oldpeak: ST depression induced by exercise relative to rest
#num: diagnosis of heart disease (angiographic disease status)
  #Value 0: < 50% diameter narrowing
  #Value 1: > 50% diameter narrowing

heart.completo <- read.csv("C:\\Users\\gabi_\\Documents\\USP\\Análise Multivariada\\Monitoria\\Aula AD\\data.csv")
heart.completo <- na_if(heart.completo,"?")
heart <- heart.completo[,-c(2,3,6,7,9,11,12,13)] #Selecionando variáveis numéricas
str(heart)
heart <- na.omit(heart) #Retirando observações faltantes
heart[,-6] <- apply(heart[,-6],2,as.numeric)
heart$num <- as.factor(heart$num) #Variável resposta

#Verificando homocedasticidade
boxM(heart[,1:5],heart[,6])
#Hipótese de homocedasticidade rejeitada

fitq <- qda(num ~ age + trestbps + chol + thalach + oldpeak, data = heart, method="mle")
fitq$scaling

fitq.escores <- predict(fitq,heart[,1:5])
ctq <- table(heart$num, fitq.escores$class)
ctq
sum(diag(prop.table(ctq))) # prop. classificação correta total

#Na prática deve-se utilizar validação cruzada
set.seed(421)
heart$id <- 1:nrow(heart)
train <- heart %>% dplyr::sample_frac(.75)
test  <- anti_join(heart, train, by = 'id')

fit.train <- qda(num ~ age + trestbps + chol + thalach + oldpeak, data = train, method="mle",prior=c(1,1)/2)
escores.test <- predict(fit.train,test[,1:5])
ct.test <- table(test$num, escores.test$class)
ct.test
sum(diag(prop.table(ct.test))) # prop. classificação correta total

#Normalidade????
#Página 595 Livro Johnson
#Páginas 177 e 186 do Livro Johnson - como verificar normalidade multivariada

dist0 <- mahalanobis(heart[heart$num == "0",c(1:5)],
                     colMeans(heart[heart$num == "0",c(1:5)]),
                     cov(heart[heart$num == "0",c(1:5)]))
dist1 <- mahalanobis(heart[heart$num == "1",c(1:5)],
                     colMeans(heart[heart$num == "1",c(1:5)]),
                     cov(heart[heart$num == "1",c(1:5)]))

#Distância de Mahalanobis
par(mfrow=c(1,2))
#df = p
qqplot(qchisq(ppoints(500), df = 5), dist0, xlab = "Quantis Teóricos", ylab = 'Distância de Mahalanobis para Grupo 0')
abline(0, 1, col = 'gray')
qqplot(qchisq(ppoints(500), df = 5), dist1, xlab = "Quantis Teóricos", ylab = 'Distância de Mahalanobis para Grupo 1')
abline(0, 1, col = 'gray')
par(mfrow=c(1,1))

#Teste
mvn(heart[heart$num == "0",c(1:5)], mvnTest = "mardia")
mvn(heart[heart$num == "1",c(1:5)], mvnTest = "mardia")

#Possíveis soluções?

##### Análise discriminante via regressão logística ######
#sex: sex 
  #1: male
  #0: female)
#cp: chest pain type
  #1: typical angina
  #2: atypical angina
  #3: non-anginal pain
  #4: asymptomatic
#fbs: (fasting blood sugar > 120 mg/dl) 
  #1: true
  #0: false)
#restecg: resting electrocardiographic results
  #0: normal
  #1: having ST-T wave abnormality (T wave inversions and/or ST elevation or depression of > 0.05 mV)
  #2: showing probable or definite left ventricular hypertrophy by Estes' criteria
#exang: exercise induced angina 
  #1: yes
  #0: no

heart2 <- heart.completo[,-c(11:13)]
str(heart2)
heart2 <- na.omit(heart2) #Retirando observações faltantes
heart2[,-c(2,3,6,7,9,11)] <- apply(heart2[,-c(2,3,6,7,9,11)],2,as.numeric)
for(i in c(2,3,6,7,9,11)) heart2[,i] <- as.factor(heart2[,i])
str(heart2)

fitrl <- glm(num ~ .,data=heart2, family=binomial(link="logit")) 
fitrl
summary(fitrl)

fitrl.logit <- predict.glm(fitrl, heart2[,-11]) #predição do logito: ln(p/(1-p))
plot(fitrl.logit, main="AD via regressão logística", pch=23, bg=c('red','green')[as.factor(fitrl.logit >= 0)])

ctrl <- table(heart2$num,fitrl.logit >= 0) #tabela com as classificações 
ctrl
sum(diag(prop.table(ctrl))) # prop. classificação correta total


########
cancer <- read.csv("C:\\Users\\gabi_\\Documents\\USP\\Análise Multivariada\\Monitoria\\Aula AD\\data_breast_cancer.csv")
cancer<- cancer[,-c(1,33)]
str(cancer)

boxM(cancer[,-1],cancer$diagnosis)

dim(cancer)
fit.cancer <- qda(diagnosis ~ ., data = cancer)
fit.cancer$scaling

escore.cancer <- predict(fit.cancer, cancer[,-1])
table(escore.cancer$class, cancer$diagnosis)

#Normalidade
distM <- mahalanobis(cancer[cancer$diagnosis == "M",-1],
                     colMeans(cancer[cancer$diagnosis == "M",-1]),
                     cov(cancer[cancer$diagnosis == "M",-1]))
distB <- mahalanobis(cancer[cancer$diagnosis == "B",-1],
                     colMeans(cancer[cancer$diagnosis == "B",-1]),
                     cov(cancer[cancer$diagnosis == "B",-1]))

#Distância de Mahalanobis
par(mfrow=c(1,2))
#df = p
qqplot(qchisq(ppoints(500), df = 30), distM, xlab = "Quantis Teóricos", ylab = 'Distância de Mahalanobis para Grupo 0')
abline(0, 1, col = 'gray')
qqplot(qchisq(ppoints(500), df = 30), distB, xlab = "Quantis Teóricos", ylab = 'Distância de Mahalanobis para Grupo 1')
abline(0, 1, col = 'gray')
par(mfrow=c(1,1))

#Teste
mvn(heart[heart$num == "0",c(1:5)], mvnTest = "mardia")
mvn(heart[heart$num == "1",c(1:5)], mvnTest = "mardia")

fit.lin.cancer <- lda(diagnosis ~. , data=cancer)
names(fit.lin.cancer)
fit.lin.cancer$scaling

escores.lin.cancer <- predict(fit.lin.cancer, cancer[,-1])
escores.lin.cancer$x

par(mfrow=c(1,2))
plot(escores.lin.cancer$x, col = cancer$diagnosis)
plot(escores.lin.cancer$x, col = escores.lin.cancer$class)

table(escores.lin.cancer$class, cancer$diagnosis)

#Com VC
fit.lin.cancer.vc <- lda(diagnosis ~. , data=cancer, CV = T)
table(cancer$diagnosis, fit.lin.cancer.vc$class)