###############################################################################
## Inhalt:
## Kapitel 7.2 (Elemente Multivariater Statistik) des Kurses von 
## Ruckdeschel & Kohl
###############################################################################

###############################################################################
## Hauptkomponenten- und Faktoranalyse
## Lesen Sie Abschnitt 7.2.4 im Skript von Ruckdeschel & Kohl
###############################################################################

###########################################################
## Hauptkomponnentenanalyse (PCA: principal component analysis)
###########################################################
## 1. Lade Datensatz
data(USArrests)

## 2. Berechne PCA 
pc11 <- princomp(USArrests) # basiert auf "eigen"
str(pc11) # Struktur des Ergebnisses

## Vergleiche
pc11$sdev^2/49*50 # umrechnen von Division durch n auf n-1
pc11$loadings[,]
eigen(cov(USArrests))

## bzw.
pc12 <- prcomp(USArrests) # basiert auf "svd" - höhere Genauigkeit
str(pc12) # Struktur des Ergebnisses

## 3. Ausgabe der Ergebnisse
summary(pc11)
loadings(pc11) # die leeren Einträge sind nicht wirklich leer bzw. 0, sondern nur klein
## bzw.
pc12
summary(pc12)

#######################################
## Sehr hoher Erklärungswert der ersten Hauptkomponente ... gut?
## Hier: Nein!
## Achtung: Hier ist die Varianz der Variablen sehr unterschiedlich!
var(USArrests)
## Folglich dominiert die Variable "Assault"
## Ausweg: Skalierung vor Anwendung der PCA

## 4. Berechne PCA mit Skalierung
pc21 <- princomp(USArrests, cor = TRUE)
pc22 <- prcomp(USArrests, scale = TRUE)

## 5. Ausgabe der Ergebnisse
summary(pc21)
## vgl.
pc21$sdev^2
pc21$loadings[,]
eigen(cor(USArrests))

loadings(pc21)
## bzw.
pc22
summary(pc22)

## 6. Graphische Darstellung der Ergebnisse
plot(pc21) # sog. screeplot
## bzw.
plot(pc22)

## oder auch
biplot(pc21)
## bzw.
biplot(pc22)

## 7. Es gibt auch ein Formelinterface im Fall beider Funktionen
princomp(~ Murder + Assault + UrbanPop + Rape, data = USArrests)
## oder kürzer
princomp(~ ., data = USArrests, cor = TRUE)

## bzw.
prcomp(~ Murder + Assault + UrbanPop + Rape, data = USArrests, scale = TRUE)
## oder kürzer
prcomp(~ ., data = USArrests, scale = TRUE)


###########################################################
## Ergänzungen mittels anderer R-Pakete
###########################################################

#######################################
## Paket "chemometrics"
#######################################
library(chemometrics)

## Kreuzvalidierung der "erklärten Varianz"
res <- pcaCV(USArrests, segments=4, repl=100, cex.lab=1.2, ylim=c(0,1), las=1)

## Erklärte Varianz durch die ersten beiden Komponenten
res <- pcaVarexpl(USArrests, a=2)

## Ausreißerdiagnostik
res <- pcaDiagplot(USArrests, pc21, a=2)

## Dies sollte aber besser mit einer robusten Kovarianzschätzung
## erfolgen!
require(robustbase)

## Robuste Kovarianzschätzung mit dem sog. Minimum Covariance Determinant
## Schätzer
(rob.cov <- covMcd(USArrests, cor = TRUE))
cov(USArrests)

## Robuste PCA
rob.pca <- princomp(USArrests, cor = TRUE, covmat = rob.cov)

## Zum Vergleich
par(mfrow = c(1, 2))
biplot(pc21, main = "Klassische Kovarianzschätzung")
biplot(rob.pca, main = "Robuste Kovarianzschätzung")

## Keine Ausreißer erkennbar
res <- pcaDiagplot(USArrests, rob.pca, a=2)

## Zum Vergleich
## ohne Skalierung ist ein Ausreißer zu erkennen
rob.pca1 <- princomp(USArrests, covmat = rob.cov)
res <- pcaDiagplot(USArrests, rob.pca1, a=2)

which.max(res$SDist)

#######################################
## Pakete "vegan" und "BiodiversityR"
#######################################
library(vegan)
library(BiodiversityR)
rda1 <- rda(USArrests)
PCAsignificance(rda1)
plot.rda <- ordiplot(rda1, choices=c(1,2), scaling=1)
ordiequilibriumcircle(rda1, plot.rda)


#######################################
## Paket "ade4"
#######################################
library(ade4)
dudi1 <- dudi.pca(USArrests, scan = FALSE, nf = 4)
dudi1$c1


#######################################
## Paket "FactoMineR"
#######################################
library(FactoMineR)
pc3 <- PCA(USArrests)
pc3
pc3$eig
pc3$var


#######################################
## Paket "pcaPP"
#######################################
## robuste PCA
library(pcaPP)
pc4 <- PCAproj(USArrests, k = 4, scale = mad)
summary(pc4)
loadings(pc4)
## zum Vergleich
summary(pc21)
loadings(pc21)

## angelehnt an princomp
## S3-Klasse "pcaPP" abgeleitet von S3-Klasse "princomp"
class(pc4)

## Graphische Darstellungen
par(mfrow = c(1, 2))
plot(pc4, main = "Robuste PCA")
plot(pc21, main = "Klassische PCA")

par(mfrow = c(1, 2))
biplot(pc4, main = "Robuste PCA")
biplot(pc21, main = "Klassische PCA")

## zur besseren Vergleichbarkeit
par(mfrow = c(1, 2))
biplot(pc4, main = "Robuste PCA", xlabs = rownames(USArrests))
biplot(pc21, main = "Klassische PCA")


#######################################
## Paket "rrcov"
#######################################
## robuste PCA
library(rrcov)
pc5 <- PcaHubert(USArrests, k = 4) # basiert auf der Kovarianz!
pc5@loadings
pc12

## Ausreißeridentifikation
plot(pc5)
## bzw.
plot(PcaHubert(USArrests, k = 2))

## Zum Vergleich
par(mfrow = c(1, 2))
biplot(pc5, main = "Robuste PCA")
stats::biplot(pc11, main = "Klassische PCA")


###############################################################################
## Übungsaufgabe 1:
## Verwenden Sie den Datensatz "banknote" aus dem Paket "dr".
## (a) Führen Sie eine (klassische) PCA durch, wobei Sie die Variable "Y" 
## bitte nicht berücksichtigen. Stellen Sie das Ergebnis geeignet graphisch dar.
## Im bi-plot verwenden Sie für die echten Banknoten (d.h., Y = 0) als Symbol
## die Null und für die gefälschten Banknoten (d.h., Y = 1) die Eins.
## Warum sollte man im vorliegenden Fall die Daten vor der PCA-Analyse skalieren?
## (b) Führen Sie eine robuste PCA durch und vergleichen Sie das Ergebnis mit
## der klassischen PCA aus Teil (a).
## (c) Verwenden Sie die robuste PCA, um Ausreißer in den Daten zu finden.
###############################################################################



###########################################################
## Faktoranalyse
###########################################################
## factanal mit varimax
?factanal
?varimax
(fact1 <- factanal(~Murder+Assault+UrbanPop+Rape, factors = 1, data = USArrests))

## siehe ?ability.cov
(ability.FA <- factanal(factors = 1, covmat=ability.cov))
update(ability.FA, factors=2)
update(ability.FA, factors=2, rotation="promax")

#######################################
## Erweiterungen und Ergänzungen finden sich in den Paketen:
## FAiR, FactoMineR, ifa, ...
#######################################



###############################################################################
## Multidimensional Scaling (MDS)
## Lesen Sie Abschnitt 7.2.5 im Skript von Ruckdeschel & Kohl
##
## Ausführliche Beschreibungen von MDS
## http://www.wiwi.uni-bielefeld.de/~frohn/Mitarbeiter/Handl/mvamdszu.pdf
##
## zusätzlich noch GGvis bzw. X/GGobi
## http://www-stat.wharton.upenn.edu/~buja/PAPERS/paper-mds-jcgs.pdf
###############################################################################
## Klassisch: metrisches MDS
## d.h., ursprüngliche Abstände und MDS-Abstände möglichst ähnlich
## evtl. muss mit "dist" zuerst die Abstandsmatrix berechnet werden!!!
loc <- cmdscale(eurodist) # eurodist enthält bereits Abstände
x <- loc[,1]
y <- -loc[,2]
plot(x, y, type="n", xlab="", ylab="", main="cmdscale(eurodist)")
text(x, y, rownames(loc), cex=0.8)

cmdsE <- cmdscale(eurodist, k=20, add = TRUE, eig = TRUE, x.ret = TRUE)
utils::str(cmdsE)

## isoMDS: nicht-metrisches MDS
## d.h., ursprüngliche Abstände und MDS-Abstände mit gleicher Anordnung
library(MASS)
loc1 <- isoMDS(eurodist)
plot(loc1$points[,1], -loc1$points[,2], type = "n")
text(loc1$points[,1], -loc1$points[,2], labels = rownames(loc1$points), cex=0.8)

## Zum Vergleich
loc
loc1

## diagnositischer Plot
## Vergleich der ursprünglichen Abstände mit den MDS-Abständen
shep <- Shepard(eurodist, loc1$points)
plot(shep, pch = 20)
lines(shep$x, shep$yf, type = "S")


## Interface zu isoMDS
library(labdsv)
loc3 <- nmds(eurodist)
plot(loc3)
plotid(loc3, ids=rownames(loc3$points)) # interaktives Identifizieren


## metaMDS basiert auf isoMDS
sol <- metaMDS(swiss) # benötigt ursprünglichen data.frame!
plot(sol, type="t")


## sammon: passt kleine Distanzen besser an
loc2 <- sammon(eurodist) # aus MASS
plot(loc2$points[,1], -loc2$points[,2], type = "n")
text(loc2$points[,1], -loc2$points[,2], labels = rownames(loc2$points), cex=0.8)

## Zum Vergleich
par(mfrow = c(1, 2))
plot(loc1$points[,1], -loc1$points[,2], type = "n", main = "isoMDS")
text(loc1$points[,1], -loc1$points[,2], labels = rownames(loc1$points), cex=0.8)
plot(loc2$points[,1], -loc2$points[,2], type = "n", main = "sammon")
text(loc2$points[,1], -loc2$points[,2], labels = rownames(loc2$points), cex=0.8)


###############################################################################
## Übungsaufgabe 2:
## Verwenden Sie den Datensatz "banknote" aus dem Paket "dr".
## (a) Lassen Sie die Variable Y aus dem Datensatz "banknote" weg und 
## berechnen Sie metrisches und nicht-metrisches MDS. 
## (b) Stellen Sie das Ergebnis graphisch dar, wobei Sie die Werte, die zu 
## echten Banknoten (d.h., Y = 0) gehören grün einfärben und Werte, die zu
## gefälschten Banknoten gehören (d.h., Y = 1) rot einfärben.
###############################################################################



###############################################################################
## Lösungen:
###############################################################################

###########################################################
## Übungsaufgabe 1
###########################################################
#######################################
## Teil (a)
#######################################
library(dr)
data(banknote)

diag(cov(banknote[,-7]))
## Die Streuung/Varianz der Variablen ist sehr unterschiedlich. Daher empfiehlt
## es sich vor der PCA eine Skalierung durchzuführen.

## klassische PCA mit Skalierung
pca.kl <- prcomp(banknote[,-7], scale = TRUE)
pca.kl
summary(pca.kl)

## oder auch
library(vegan)
library(BiodiversityR)
rda.kl <- rda(banknote[,-7], scale = TRUE)
PCAsignificance(rda.kl)

## scree- und bi-plot
plot(pca.kl)
stats::biplot(pca.kl, xlabs = as.character(banknote$Y))


#######################################
## Teil (b)
#######################################
library(robustbase)
rob.cov <- covMcd(banknote[,-7], cor = TRUE)

## Robuste PCA
pca.rob <- princomp(banknote[,-7], cor = TRUE, covmat = rob.cov)

## Zum Vergleich
par(mfrow = c(1, 2))
stats::biplot(pca.kl, main = "Klassische Kovarianzschätzung", xlabs = as.character(banknote$Y))
stats::biplot(pca.rob, main = "Robuste Kovarianzschätzung", xlabs = as.character(banknote$Y))


#######################################
## Teil (c)
#######################################
## Ausreißer?
library(chemometrics)
res <- pcaDiagplot(banknote[,-7], pca.rob, a=2, pch = as.character(banknote$Y))

## bzw. (basiert auf Kovarianz, d.h., keine Skalierung!)
library(rrcov)
plot(PcaHubert(banknote[,-7], k = 6))


###########################################################
## Übungsaufgabe 2
###########################################################

#######################################
## Teil (a)
#######################################
library(dr)
data(banknote)

banknote.dist <- dist(banknote[,-7])
## metrisch
mds.kl <- cmdscale(banknote.dist)

## nicht-metrisch
mds.iso <- isoMDS(banknote.dist)
mds.sam <- sammon(banknote.dist)

#######################################
## Teil (b)
#######################################
plot(mds.kl, type="p", pch = 20, xlab="", ylab="", 
     col = ifelse(banknote$Y, "red", "green"), main="cmdscale")
plot(mds.iso$points, type="p", pch = 20, xlab="", ylab="", 
     col = ifelse(banknote$Y, "red", "green"), main="isoMDS")
plot(mds.sam$points, type="p", pch = 20, xlab="", ylab="", 
     col = ifelse(banknote$Y, "red", "green"), main="sammon")