#  ------------------------------------------------------------------------
# Diskriminacni analyza ---------------------------------------------------
#  ------------------------------------------------------------------------

install.packages("ade4")
install.packages("klaR")
install.packages("ellipse")
install.packages("ggplot2")
install.packages("scales")
install.packages("gridExtra")

library("ade4")
library("klaR")
library("ellipse")
library("ggplot2")
library("scales")
library("gridExtra")

setwd("c:/Users/brozoval/Desktop/6. cvičení - DA/data/")
getwd()

############################################################################################################
#############################################    IRIS dataset       ########################################
############################################################################################################

data(iris)# nacte data kosatcu (sepal = kalisni listek, petal = korunni listek)
head(iris, 3)
names(iris)
dim(iris)
table(iris$Species)

par(mar = c(0, 0, 0, 0))
pan1 <- function(x, y, ...) {
xy <- cbind.data.frame(x, y)
s.class(xy, iris$Species, include.ori = F, add.p = T, clab = 1.5,
col = c("blue", "black", "red"), cpoi = 2, csta = 0.5)
}
pairs(iris[, 1:4], panel = pan1)
# vidime, ze delka a sirka korunnich listku nejlepe diskriminuje skupinu setosa

# vyber promennych pomoci dopredne selekce 
?greedy.wilks
iris.step<- greedy.wilks(Species ~  Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, data=iris, niveau = 0.1)
iris.step 						                  # jako prvni byla pridana promenna Petal.Length, jako 2. Sepal.Width, jako 3. Petal.Width, jako posledni Sepal.Length
							                          # vsechny 4 promenne byly vybraby do modelu, pokud by promenna nebyla vybrana, nebyla by uvedena ve vyctu
iris.step$results$Wilks.lambda  	      # Wilksovo lambda modelu - snizuje se s poctem promennych v modelu - lepsi diskriminace
iris.step$results$F.statistics.overall	# F statistika pro cely model s aktualne vybranymi promennymi 
iris.step$results$p.value.overall		    # jeji statisticka vyznamnost
iris.step$results$F.statistics.diff		  # F statistika - srovnava model s a bez nove promenne 
iris.step$results$p.value.diff		      # jeji statisticka vyznamnost (p < 0.05 -> promenna vyznamne meni diskriminaci modelu)

# dopredna selekce vybrala do modelu vsechny 4 promenne, finalni F statistika modelu = 199 s p = 1.365006e-112 a celkovym Wilksovym lambda modelu = 0.0234

# vypocet diskriminacni analyzy pro vybrane promenne
?lda
DA <- lda(formula = Species ~ .,data = iris,prior = c(1,1,1)/3) 
DA
				                # . rika vezmi vsechny ostatni parametry ze souboru iris jako prediktory
DA$prior		          	# apriorni pravedpodobnost - rovnomerna
DA$counts			          # pocty objektu ve skupinach
DA$means			          # prumery v ramci skupin (muzeme vycist, ktere skupiny se lisi v jakem parametru, napr. setosa ma nejkratsi okvetni listky)
DA$scaling		        	# vlastni vektory - linearni kombinace puvodnich promennych, dle kterych je vytvorena nova osa (vzdy maximalne pocet skupin/promennych minus jedna), 

(DA$svd^2/sum(DA$svd^2))# podil mezi a vnitro shlukove variability na osach; prvni osa popisuje vice nez 99 % mezishlukove variability dat

# Srovnani skutecneho typu kosatce a klasifikace pomoci diskriminacni analyzy
names(predict(DA))                # $class = klasifikace objektu do skupin, $posterior = posteriorni pst, $x = pozice objektu na novych osach
iris_table<-table(Skutecnost=(Skutecnost<-iris$Species),Klasifikace=(Klasifikace<-predict(DA)$class))
iris_table					              # 1 virginica byla chybne klasifikovana jako versicolor a 2 versicilor byly klasifikovany jako virginica
diag(iris_table)			 	          # pocty spravne zaklasifikovanych  pripadu pro jednotlive skupiny
diag(prop.table(iris_table,1))	  # procento spravne zaklasifikovanych  pripadu pro jednotlive skupiny
sum(diag(iris_table))/nrow(iris)	# presnost klasifikace

# Matice posteriornich pravdepodobnosti se zarazenim objektu do skupin
post.pst<-predict(DA)$posterior		# posteriorni pravdepodobnost pro objekt, ze patri do dane skupiny 
klas<-predict(DA)$class						# klasifikace objektu
result<-as.data.frame(cbind(post.pst,as.character(klas),as.character(iris[,5])))	# matice posteriornich pravdepodobnosti s klasifikaci do skupin
colnames(result)[4:5]<-c("DA","puvodni")
result
result[which(result$DA!=result$puvodni),]# u techto objektu byla chybne provedena klasifikace

# Graf pro objekty ve vicerozmernem prostoru kanonickych os obarvenych podle prislusnosti objektu do skupiny

pozice<-predict(DA)$x				# pozice objektu na novych osach

windows()
par(mfrow=c(1,2))
plot(pozice[,1], pozice[,2], type="n",main="nova klasifikace")
text(pozice[,1], pozice[,2], rownames(iris), col=c("blue", "black", "red")[as.numeric(klas)])   # "nova" klasifikace
legend(0,3,col=c("blue","black","red"),legend=c("setosa","versicolor","virginica"),lty=c(1,1,1))
abline(h=0, lty="dotted")
abline(v=0, lty="dotted")

# vykresleni elips kolem skupin bodu 
for (i in 1:length(levels(klas))){
  cov <- cov(pozice[as.numeric(klas)==i,])
  centre <- apply (pozice[as.numeric(klas)==i,],2,mean)
  lines(ellipse(cov, centre=centre, level=0.95))
}

plot(pozice[,1], pozice[,2], type="n",main="puvodni klasifikace")
text(pozice[,1], pozice[,2], rownames(iris), col=c("blue", "black", "red")[as.numeric(iris[,5])])   # puvodni klasifikace
abline(h=0, lty="dotted")
abline(v=0, lty="dotted")
legend(0,3,col=c("blue","black","red"),legend=c("setosa","versicolor","virginica"),lty=c(1,1,1))

# vykresleni elips kolem skupin bodu 
for (i in 1:length(levels(iris[,5]))){
  cov <- cov(pozice[as.numeric(iris[,5])==i,])
  centre <- apply (pozice[as.numeric(iris[,5])==i,],2,mean)
  lines(ellipse(cov, centre=centre, level=0.95))
}

# 71, 84 klasifikovany jako virginica (vlevo v grafu cervene), ale ve skutecnosti patri do versicolor (vpravo cerne)
# 134 klasifikovan do versicolor (vlevo v grafu cerne), ale ve skutecnosti patri do virginica (vpravo cervene)

############################################################################################################
#########################################          validace          #######################################
############################################################################################################

#krosvalidace 1:1

train <- sample(1:150, 75)							# nahodne vybereme trenovaci soubor jako polovinu pozorovani z datasetu, 
											                  # tento prikaz vytvori indexy, podle kterych budou pozorovani vybrana
											                  # pozor! Pri opetovnem spusteni budou vygenerovany nove indexy a bude davat odlisne vysledky
table(iris[train,5])								    # vidime, ze i na podsouboru jsou skupiny rovnomerne zastoupeny
DA3 <- lda(Species ~ .,iris, prior = c(1,1,1)/3,subset = train)	# tvorba modelu na trenovacim souboru
 
plda = predict(object = DA3,newdata = iris[-train, ]) 		# klasifikace na zbytku
plda$class 

# Srovnani skutecneho typu kosatce a klasifikace pomoci diskriminacni analyzy na testovaci souboru (ten, ktery nebyl vyuzit pro tvorbu modelu)
iris_table<-table(Skutecnost=(Skutecnost<-iris$Species[-train]), 
      Klasifikace=(Klasifikace<-plda$class))
iris_table
diag(iris_table)			                  # pocty spravne zaklasifikovanych  pripadu pro jednotlive skupiny
diag(prop.table(iris_table,1))          # procento spravne zaklasifikovanych  pripadu pro jednotlive skupiny
sum(diag(iris_table))/nrow(iris[-train,])# presnost klasifikace

############################################################################################################
#########################################    srovnani PCA a DA       #######################################
############################################################################################################

# DA natoci prostor tak, aby nejvice odlisila pozice subjektu mezi skupinami
# PCA natoci prostor tak, aby popsala maximum puvodniho rozptylu dat

pca <- prcomp(iris[,-5],center = TRUE,scale. = TRUE) 
prop.pca = pca$sdev^2/sum(pca$sdev^2)

lda <- lda(Species ~ .,  iris, prior = c(1,1,1)/3)
prop.lda = DA$svd^2/sum(DA$svd^2)
plda <- predict(object = lda,newdata = iris)

dataset = data.frame(species = iris[,"Species"],pca = pca$x, lda = plda$x)

p1 <- ggplot(dataset) + geom_point(aes(lda.LD1, lda.LD2, colour = species, shape = species), size = 2.5) + 
  labs(x = paste("LD1 (", percent(prop.lda[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda[2]), ")", sep="")) 

p2 <- ggplot(dataset) + geom_point(aes(pca.PC1, pca.PC2, colour = species, shape = species), size = 2.5) +
  labs(x = paste("PC1 (", percent(prop.pca[1]), ")", sep=""),
       y = paste("PC2 (", percent(prop.pca[2]), ")", sep=""))

grid.arrange(p1, p2)

#################################################################################################
#################################        TITANIC     ############################################
#################################################################################################

dataDA.load<-read.table('titanic.csv', sep=";", header=TRUE, dec=".")

dim(dataDA.load)
names(dataDA.load)
table(dataDA.load$preziti)
head(dataDA.load)
tail(dataDA.load)

dataDA<-dataDA.load[1:1040,]
dim(dataDA)

# vztah predktoru k preziti pomoci logisticke regrese
summary(glm(preziti~pohlavi,family = "binomial",data=dataDA))
summary(glm(preziti~vek,family = "binomial",data=dataDA))
summary(glm(preziti~pocet.sourozencu.partneru,family = "binomial",data=dataDA))
summary(glm(preziti~pocet.deti.rodicu,family = "binomial",data=dataDA))
summary(glm(preziti~cenalistku,family = "binomial",data=dataDA))
summary(glm(preziti~trida,family = "binomial",data=dataDA))

### vypocet diskriminacni analyzy
# vyber parametru 
DA.step<- greedy.wilks(preziti~
                         pohlavi+
                         vek+
                         pocet.sourozencu.partneru+
                         pocet.deti.rodicu+
                         cenalistku+
                         trida,data=dataDA, niveau = 0.05)
DA.step

# nastaveni apriornich pravdepodobnosti
preziti.tab<-table(dataDA[,"preziti"])
preziti.tab
apriori<-c(preziti.tab/1040)
apriori

# model DA
DA <- lda(formula = preziti~
            pohlavi+
            trida+
            vek+
            pocet.sourozencu.partneru+
            cenalistku,
          data = dataDA,
          prior = apriori) 
DA

# skutecnost vs. predikce
class_table<-table(Skutecnost=(Skutecnost<-dataDA$preziti), 
                   Klasifikace=(Klasifikace<-predict(DA)$class))
class_table			                        
diag(class_table) 	                    
1040-sum(diag(class_table))
diag(prop.table(class_table,1))
sum(diag(class_table))/nrow(dataDA)

# matice posteriornich pravdepodobnosti s klasifikaci do skupin
post.pst<-predict(DA)$posterior
klas<-predict(DA)$class
result<-as.data.frame(cbind(post.pst,as.character(klas),as.character(dataDA[,"preziti"])))
colnames(result)[3:4]<-c("DA","puvodni")
head(result)
result[which(result$DA!=result$puvodni),]
result[c(932,1009),]

# klasifikace na zbytku
plda = predict(object = DA,newdata = dataDA.load[1041:1045, ])	
table(plda$class)

#################################################################################################
#################################       Znalosti     ############################################
#################################################################################################

dataDA.load<-read.table('znalosti.txt', sep="\t", header=TRUE, dec=".")
dim(dataDA.load)
dataDA.load$znalosti
names(dataDA.load)

dataDA<-dataDA.load[1:235,]

head(dataDA)
dim(dataDA)
str(dataDA)
summary(dataDA)
names(dataDA)

# vztah prediktoru k endpointu
windows()
par(mfrow=c(3,3))
boxplot(dataDA$Biologie~as.factor(dataDA$znalosti),main="znalosti~ Biologie",ylab="Biologie",xlab="znalosti")
boxplot(dataDA$matematika~as.factor(dataDA$znalosti),main="znalosti~ matematika",ylab="matematika",xlab="znalosti")
boxplot(dataDA$cestina~as.factor(dataDA$znalosti),main="znalosti~ cestina",ylab="cestina",xlab="znalosti")
boxplot(dataDA$cizijazyk~as.factor(dataDA$znalosti),main="znalosti~ cizijazyk",ylab="cizijazyk",xlab="znalosti")
boxplot(dataDA$chemie~as.factor(dataDA$znalosti),main="znalosti~ chemie",ylab="chemie",xlab="znalosti")
boxplot(dataDA$ZSV~as.factor(dataDA$znalosti),main="znalosti~ ZSV",ylab="ZSV",xlab="znalosti")
boxplot(dataDA$fyzika~as.factor(dataDA$znalosti),main="znalosti~ fyzika",ylab="fyzika",xlab="znalosti")
boxplot(dataDA$telocvik~as.factor(dataDA$znalosti),main="znalosti~ telocvik",ylab="telocvik",xlab="znalosti")
boxplot(dataDA$dejepis~as.factor(dataDA$znalosti),main="znalosti~ dejepis",ylab="dejepis",xlab="znalosti")

### vypocet diskriminacni analyzy
# vyber parametru 
DA.step<- greedy.wilks(znalosti~., data=dataDA, niveau = 0.05)
DA.step

# nastaveni apriornich pravdepodobnosti
znalosti.tab<-table(dataDA[,"znalosti"])
apriori<-c(znalosti.tab/235)
apriori

# model DA
DA <- lda(formula = znalosti~.,data = dataDA,prior = apriori) 
DA

# skutecnost vs. predikce
class_table<-table(Skutecnost=(Skutecnost<-dataDA$znalosti), 
                   Klasifikace=(Klasifikace<-predict(DA)$class))
class_table 
diag(class_table)
235-sum(diag(class_table))
diag(prop.table(class_table,1))
sum(diag(class_table))/nrow(dataDA)

#kanonicke osy
DA$scaling		        
(DA$svd^2/sum(DA$svd^2))

# objekty na prvnich dvou osach
lda <- lda(znalosti ~ .,  dataDA, prior = apriori)
prop.lda = DA$svd^2/sum(DA$svd^2)
plda <- predict(object = lda,newdata = dataDA)
dataset = data.frame(species = as.factor(dataDA[,"znalosti"]), lda = plda$x)

p1 <- ggplot(dataset) + geom_point(aes(lda.LD1, lda.LD2, colour = species,shape = species), size = 2.5) + 
  labs(x = paste("LD1 (", percent(prop.lda[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda[2]), ")", sep="")) 
p1

# matice posteriornich pravdepodobnosti s klasifikaci do skupin
post.pst<-predict(DA)$posterior
klas<-predict(DA)$class
result<-as.data.frame(cbind(post.pst,as.character(klas),as.character(dataDA[,"znalosti"])))
colnames(result)[6:7]<-c("DA","puvodni")
head(result)
result[which(result$DA!=result$puvodni),]

# klasifikace na zbytku
plda = predict(object = DA,newdata = dataDA.load[236:239, ])
table(plda$class)