# Stochastik IV, Übungsblatt 5

# Aufgabe 3

library(lattice)

# Einlesen der Daten und Einbinden in Suchpfad

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

attach(data.olives)

# Hauptkomponentenanalyse (basierend auf der Korrelationsmatrix)

olives <- data.olives[, 3:10]  # die acht Fettsäuren

olives.pca <- princomp(olives, cor = TRUE)  

# Varianzkomponenten

names(olives.pca$sdev) <- c(paste("PC", 1:8, sep = ""))
summary(olives.pca)

# Screeplot

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

windows()

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

# Kriterium i, Eigenwertkriterium (siehe Aufgabe 4, (a))

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 = olives.pca)

# Kriterium ii (siehe Aufgabe 4, (a))

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 = olives.pca, c = 0.90)
crit.ii(name.pca = olives.pca, c = 0.95)
crit.ii(name.pca = olives.pca, c = 0.975)
crit.ii(name.pca = olives.pca, c = 0.99)
crit.ii(name.pca = olives.pca, c = 0.9975)
crit.ii(name.pca = olives.pca, c = 0.999)
crit.ii(name.pca = olives.pca, c = 0.9999)

# Kriterium iii, Scree-Test / Ellenbogen-Kriterium (siehe Aufgabe 4, (a)): 
# Inspektion von Screeplot liefert Knickstellen bei Hauptkomponente 3 oder 4, 
# mit zugehörigen Eigenwerten 1.016355435 bzw. 0.792898832. 
# Je nachdem, ergibt sich also k = 3 oder 4.

# Ladungen und Scatterplot-Matrix der Ladungen in 2D Eigenvektor-Ebenen

dimnames(olives.pca$loadings) <- list(names(olives),paste("PC", 1:8, sep = ""))
loadings(olives.pca)  # oder: olives.pca$loadings

windows()

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

# Scatterplots für die ersten drei Hauptkomponenten

windows()

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

pairs(olives.pca$scores[, 1:3], main = "SPLOM der ersten 3 Hauptkomponenten (Farbe nach Region)",
      col = as.numeric(Region), cex = 1.5, pch = 19)
title(sub = "North: black; Sardinia: red; South: green", cex.sub = .75) 

windows()

splom(~ olives.pca$scores[, 1:3] | Region, main = "Trellis SPLOM-Plot", xlab = "", cex = 0.5,
      col = "lightblue", pscales = 0, layout = c(3,1))