# Wstepna analiza danych wielowym ~ Interesujące rzuty ~ Redukcja wymiaru.
# Program wyznacza 3 okienka zawierające po 6 rysunków. Okienko 1 pokazuje
# podział danych na rzeczywiste klasy na płaszczyźnie wyznaczonej przez PCA
# i skalowanie Sammona dla oryginalnych danych, danych znormalizowanych
# oraz danych przekształconych logarytmicznie. Okienko 2 pokazuje to samo dla
# podziału wyznaczonego przez k-means, a  Okienko 3 - przez hclust.
# Możemy obserwować efekty skalowania danych, róznice w rzutowaniu PCA
# i Sammona oraz rezultaty klasteryzacji.
# Dla danych iris - analiza jest udana, dla crabs - nieudana.

# Wprowadzenie danych
#library(MASS); library(nlme)
#C=crabs[,1:2]; Kl=ifelse(C[,1]=="O",2,1) #Kl=ifelse(C[,2]=="M",2,1)
#X=as.matrix(crabs[,4:8]); nrKl=2
Kl=c(rep(1,50), rep(2,50), rep(3,50)); X=as.matrix(iris[-143,-5]);  nrKl=3
# Sprawdzimy wpływ skalowania danych na wyniki analizy
X2=t(t(X)/apply(X,2,sd)); X3=log(X)

# Analiza Składowych Głównych (Principal Component Analysis): ortogonalizujemy
# kolumny X (cechy) i ustawiamy je w porzadku malejących wariancji.
pca=prcomp(X); pca2=prcomp(X2); pca3=prcomp(X3)

# Skalowanie wielowymiarowe (Multidimensional Scaling): rzutujemy obserwacje
# na płaszczyznę (podprzestrzeń) najlepiej zachowującą odległości między
# obserwacjami. Pomijamy przypadek minimalizacji odległości euklidesowej
# (classical multidimensional scaling = metrical scaling  = principal coordinate
# analysis), ponieważ jest równoważny PCA: cmd=cmdscale(dist(X)).
# Skalowanie Sammona minimalizuje względne odległości między danymi.
sa=sammon(dist(X),trace=FALSE)$points
sa2=sammon(dist(X2),trace=FALSE)$points
sa3=sammon(dist(X3),trace=FALSE)$points

# Wykonanie 6 obrazków na 1 okienku wymaga zmiany parametrów rysowania.
oldpar=par(no.readonly=TRUE)

par(las=1, mfcol=c(3,2), oma=c(0,0,1,0))
klasy=Kl
plot(pca$x[,1], pca$x[,2], main="original data", col=klasy)
plot(pca2$x[,1],pca2$x[,2],main="scaled data",   col=klasy)
plot(pca3$x[,1],pca3$x[,2],main="log data",      col=klasy)
plot(sa[,1], sa[,2], main="original data", col=klasy)
plot(sa2[,1],sa2[,2],main="scaled data",   col=klasy)
plot(sa3[,1],sa3[,2],main="log data",      col=klasy)
title("True classes",outer=TRUE)

# Podział dany przez klasteryzację metodą k-średnich
klasy= kmeans(X, nrKl)$cluster
klasy2=kmeans(X2,nrKl)$cluster
klasy3=kmeans(X3,nrKl)$cluster;
x11(); par(las=1, mfcol=c(3,2), oma=c(0,0,1,0))
plot(pca$x[,1], pca$x[,2], main="original data", col=klasy)
plot(pca2$x[,1],pca2$x[,2],main="scaled data",   col=klasy2)
plot(pca3$x[,1],pca3$x[,2],main="log data",      col=klasy3)
plot(sa[,1], sa[,2], main="original data", col=klasy)
plot(sa2[,1],sa2[,2],main="scaled data",   col=klasy2)
plot(sa3[,1],sa3[,2],main="log data",      col=klasy3)
title("K-means",outer=TRUE)

# Podział dany przez aglomeracyjną klasteryzację hierarchiczną
hcl= hclust(dist(X));   klasy = as.vector(cutree(hcl,nrKl))
hcl2=hclust(dist(X2)); klasy2 = as.vector(cutree(hcl2,nrKl))
hcl3=hclust(dist(X3)); klasy3 = as.vector(cutree(hcl3,nrKl))
x11(); par(las=1, mfcol=c(3,2), oma=c(0,0,1,0))
plot(pca$x[,1], pca$x[,2], main="original data", col=klasy)
plot(pca2$x[,1],pca2$x[,2],main="scaled data",   col=klasy2)
plot(pca3$x[,1],pca3$x[,2],main="log data",      col=klasy3)
plot(sa[,1], sa[,2], main="original data", col=klasy)
plot(sa2[,1],sa2[,2],main="scaled data",   col=klasy2)
plot(sa3[,1],sa3[,2],main="log data",      col=klasy3)
title("Hierarch clust",outer=TRUE)

par(oldpar)

#x11(); plot(hcl) # dendrogram