# Stochastik IV, Übungsblatt 5

# Aufgabe 4

library(lattice)

# Einlesen der Daten und Einbinden in Suchpfad

data.crabs <- 
read.table("D:/Univ Augsburg/Lehre/SS07/VO_MultiVar/Datensaetze/crabs.txt",
head = TRUE, sep = "\t", quote = "")

attach(data.crabs)

# Hauptkomponentenanalyse (basierend auf der Kovarianzmatrix)

crabs <- data.crabs[, 4:8]  # die fünf morphologischen Messvariablen

crabs.pca <- princomp(crabs)  # ACHTUNG: (1) Defaulteinstellung für princomp() 
                              # benutzt 1/n statt 1/(n-1) zur Berechnung von 
                              # Kovarianzmatrix; 
                              # (2) In princomp() werden Datenmatrix-Spalten
                              # auf 0 zentriert

# Teilaufgabe (a)

summary(crabs.pca)

plot(crabs.pca, main = "Screeplot (\"barplot\")", 
     xlab = "Principal Components", col = "lightblue")

windows()

plot(crabs.pca, main = "Screeplot (\"lines\")", type = "lines", 
     col = "lightblue")
mtext(side = 1, line = 3, "Principal Components")

# Kriterium i, Eigenwertkriterium

crit.i <- function(name.pca) { 
	p <- length(name.pca$sdev)
	for ( r in c(1:p) ) 
   		if ( name.pca$sdev[r]^2 >= (sum(name.pca$sdev^2) / p) )
  			i <- r
    	i
}

crit.i(name.pca = crabs.pca)

# Kriterium ii

crit.ii <- function(name.pca, c) {
	p <- length(name.pca$sdev)
	for ( r in c(1:p) )
		if ( sum(name.pca$sdev[1:r]^2) >= (c * sum(name.pca$sdev^2)) )
			break
	r	
}

crit.ii(name.pca = crabs.pca, c = 0.90)
crit.ii(name.pca = crabs.pca, c = 0.95)
crit.ii(name.pca = crabs.pca, c = 0.975)
crit.ii(name.pca = crabs.pca, c = 0.99)
crit.ii(name.pca = crabs.pca, c = 0.9975)
crit.ii(name.pca = crabs.pca, c = 0.999)
crit.ii(name.pca = crabs.pca, c = 0.9999)

# Kriterium iii, Scree-Test / Ellenbogen-Kriterium: Inspektion von Screeplot 
# liefert Knickstelle bei Hauptkomponente 2, mit zugehörigem Eigenwert 
# 1.29035257. Je nachdem, entweder mit oder ohne Knickstelle, ergibt sich 
# k = 1 oder 2.

# Teilaufgabe (b)

loadings(crabs.pca)  # oder: crabs.pca$loadings

windows()

pairs(loadings(crabs.pca), main = "SPLOM der 5 Variablen in den
      2D Eigenvektor-Ebenen")

predict(crabs.pca)[1:3, ]  # oder: crabs.pca$scores[1:3, ]

# Teilaufgabe (c)

dimnames(crabs.pca$scores) <- list(NULL, paste("PC", 1:5, sep = ""))

windows()

pairs(crabs.pca$scores[, 1:3], main = "SPLOM der ersten 3 Hauptkomponenten (Farbe und Form nach den 4 Gruppen)",
      cex = 1.5, col = (sp - 1) + (2 * (sex - 1)) + 1, pch = (sp - 1) + (2 * (sex - 1)) + 1)
title(sub = "Blue Male: black; Blue Female: green; Orange Male: red; Orange Female: blue", 
      cex.sub = .75) 

# NOTIZ: "Blue Male": col=1, pch=1; "Blue Female": col=3, pch=3; 
#        "Orange Male": col=2, pch=2; "Orange Female": col=4, pch=4

sp <- as.factor(sp); levels(sp) <- c("Blue", "Orange")
sex <- as.factor(sex); levels(sex) <- c("Male", "Female")

windows()

splom(~ crabs.pca$scores[, 1:3] | sp * sex, main = "Trellis SPLOM-Plot", xlab = "", cex = 0.5,
      col = "lightblue", pscales = 0)