#  ------------------------------------------------------------------------
# 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/4. cviceni/")
getwd()

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

data(iris)# nacte data (sepal = kalisni listek, petal = okvetni 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 okvetnich 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 ktere 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
head(predict(DA))                 # $class = klasifikace objektu do skupin, $posterior = posteriorni pst
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

# 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
# v grafech: modrou setosa, cernou versicolor, cervenou virginica

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

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

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

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="")) # DA natoci prostor tak, aby nejvice odlisil pozice subjektu mezi skupinami

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)

############################################################################################################
#############################################    samostatny ukol    ########################################
############################################################################################################

# Nactete do R datovy soubor bankloan.csv. Datovy soubor popisuje 850 domacnosti, kde u 700 je uvedeno, zda si domacnost
# v minulosti sjednala pujcku. Na zaklade vybranych prediktoru (promenna vek az jine dluhy) sestavte pomoci diskriminacni analyzy
# model, ktery bude klasifikovat domacnosti do dvou kategorii: s/bez pujcky. Model postavte na prvnich 700 domacnostech a pak aplikujte na zbylych 150.
# Provedte:

# A) Vytvorte model na prvnich 700 domacnostech! Je potreba oriznout datovy soubor.

domacnosti<-read.table('bankloan.csv', sep=";", header=TRUE, dec=".")
domacnostiDA<-domacnosti[1:700,]
dim(domacnostiDA)
names(domacnostiDA)
head(domacnostiDA)
summary(domacnostiDA)

#   1) Vyhodnotte vztah prediktoru ke grupovaci promenne (pujcka v minulosti). Napr. vykreslete krabicovy graf: boxplot (prediktor~group).

windows()
par(mfrow=c(3,3))
boxplot(domacnostiDA$Vek ~ domacnostiDA$pujcka.v.minulosti,main="vek")
boxplot(domacnostiDA$vzdelani ~ domacnostiDA$pujcka.v.minulosti,main="vzdelani")
boxplot(domacnostiDA$doba.u.jednoho.zamestnavatele ~ domacnostiDA$pujcka.v.minulosti,main="doba u jednoho zamestnavatele")
boxplot(domacnostiDA$doba.na.stejnem.bydlisti ~ domacnostiDA$pujcka.v.minulosti,main="doba na stejnem bydlisti")
boxplot(domacnostiDA$prijem.domacnosti ~ domacnostiDA$pujcka.v.minulosti,main="prijem domacnosti")
boxplot(domacnostiDA$dluh.kreditni.karta~ domacnostiDA$pujcka.v.minulosti,main="dluh kreditni karta")
boxplot(domacnostiDA$jine.dluhy ~ domacnostiDA$pujcka.v.minulosti,main="jine dluhy")

aggregate(domacnostiDA$Vek ~ domacnostiDA$pujcka.v.minulosti, FUN = function(x) c(M=mean(x),Med=median(x), SD=sd(x)))
aggregate(domacnostiDA$vzdelani ~ domacnostiDA$pujcka.v.minulosti, FUN = function(x) c(M=mean(x),Med=median(x), SD=sd(x)))
aggregate(domacnostiDA$doba.u.jednoho.zamestnavatele~ domacnostiDA$pujcka.v.minulosti, FUN = function(x) c(M=mean(x),Med=median(x), SD=sd(x)))
aggregate(domacnostiDA$doba.na.stejnem.bydlisti~ domacnostiDA$pujcka.v.minulosti, FUN = function(x) c(M=mean(x),Med=median(x), SD=sd(x)))
aggregate(domacnostiDA$prijem.domacnosti ~ domacnostiDA$pujcka.v.minulosti, FUN = function(x) c(M=mean(x),Med=median(x), SD=sd(x)))
aggregate(domacnostiDA$dluh.kreditni.karta~ domacnostiDA$pujcka.v.minulosti, FUN = function(x) c(M=mean(x),Med=median(x), SD=sd(x)))
aggregate(domacnostiDA$jine.dluhy~ domacnostiDA$pujcka.v.minulosti, FUN = function(x) c(M=mean(x),Med=median(x), SD=sd(x)))

# pacienti s pujckou jsou mladsi, maji vyssi vzdelani, kratsi dobu u 
# jednoho zamestnavatele a na trvalem bydlisti, nizsi prijem a vyssi dluhy

#   2) Jake promenne budou vstupovat do modelu? Vyhodnotte korelace prediktoru navzajem (vek az jine dluhy) pomoci funkce cor. Zvazte vysledek dopredne selekce.

domacnostiDA<-domacnostiDA[,-1]         # ID domacnosti neni prediktor, radeji smazeme nebo ulozime do row.names
korelace<-cor(domacnostiDA)
korelace[lower.tri(korelace)]
max(abs(korelace[lower.tri(korelace)]))	# maximalni korelace 0.63 - jine dluhy a dluhy na kreditni karte jsou nejvice korelovane promenne v souboru
							                          # neni extremne vysoka hodnota - tedy v datech nemame redundantni promenne

names(domacnostiDA)			                # mame 7 prediktoru a jednu klasifikacni promennou: pujcka.v.minulosti
?greedy.wilks
domacnostiDA.step<- greedy.wilks(pujcka.v.minulosti~ ., data=domacnostiDA, niveau = 0.05)
domacnostiDA.step				                # vidime, ze celkove model spatne klasifikuje domacnosti do skupin s a bez pujcky (Wilksovo lambda ma vysokou hodnotu)
						                            # vzhledem k faktu, ze u zbylych domacnosti mame vyhodnoceny vsechny promenne a nechceme model aplikovat v praxi
						                            # vybereme do analyzy vsechny promenne, ktere nam dopredna selekce doporucila

#   3) Provedte diskriminacni analyzu pro vybrane promenne (tedy prvnich 700 domacnosti=radku a konkretni sloupce=vybrane prediktory ze vstupni tabulky).
#      Apriorni pravdepodobnosti nastavte proporcionalni nasemu souboru. Uvedte Wilksovo lambda pro cely model a interpretujte.

domacnostiDA.step			# Wilksovo lambda = 0.74, vysoka hodnota, tedy model spatne diskriminuje skupiny,
					            # asociovana F statistika vychazi vyznamne, muze byt ale dano velkym N

pujcka.n<-table(domacnostiDA[,ncol(domacnostiDA)])
pujcka.n
apriori<-c(pujcka.n[1]/700,pujcka.n[2]/700)
apriori								# apriorni pravdepodobnosti nastaveny proporcionalne nasemu souboru

DA <- lda(formula = pujcka.v.minulosti~ 
	doba.u.jednoho.zamestnavatele +
	dluh.kreditni.karta +
	doba.na.stejnem.bydlisti +
	jine.dluhy +
	prijem.domacnosti,
	data = domacnostiDA,prior = apriori) 
DA

#   4) Kolik dostavate kanonickych os. Uvedte v procentech popis variability mezi skupinami na novych kanonickych osach.

DA$scaling 			        # mame jednu kanonickou osu, protoze klasifikujeme 2 skupiny - 1 osa
DA$svd^2/sum(DA$svd^2)	# jelikoz mame pouze jednu osu, tato musi popisovat 100 % variability

#   5) Vyhodnotte spravnost modelu - pocet a procento spravne klasifikovanych objektu. Vypiste domacnosti (s posteriornimi pravdepodobnostmi), ktere jsou klasifikovany spatne. 

domacnostiDA_table<-table(Skutecnost=(Skutecnost<-domacnostiDA$pujcka.v.minulosti),Klasifikace=(Klasifikace<-predict(DA)$class))
domacnostiDA_table					            # 121 domacnosti s pujckou klasifikovano jako bez pujcky, 17 bez pujcky klasifikovano jako s pujckou
diag(domacnostiDA_table)			 	        # pocty spravne zaklasifikovanych  pripadu pro jednotlive skupiny
diag(prop.table(domacnostiDA_table,1))	# procento spravne zaklasifikovanych  pripadu pro jednotlive skupiny
								                        # model spatne rozezna domacnosti s pujckou, domacnosti bez pujcky klasifikuje s pomerne velkou uspesnosti

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,klas,domacnostiDA$pujcka.v.minulosti))		# matice posteriornich pravdepodobnosti s klasifikaci do skupin
colnames(result)[3:4]<-c("DA","puvodni")
head(result)
result$DA<-result$DA-1					                # R klasifikovalo objekty do skupina 1-2 namisto 0-1.
result[which(result$DA!=result$puvodni),]		    # chybne klasifikovane domacnosti
dim(result[which(result$DA!=result$puvodni),])	# 138 chybne klasifikovanych domacnosti

# B) Na poslednich 150 domacnostech aplikujte vytvoreny model a uvedte pocet domacnosti, ktere model klasifikuje bez a s pujckou.

plda = predict(object = DA,newdata = domacnosti[701:850, ])	# klasifikace na zbytku
table(plda$class)								                            # model klasifikoval 132 domacnosti bez pujcky, 18 s pujckou