library(MASS)
library(lattice)
help(fgl)
head(fgl)

n <- dim(fgl)[1]
p <- dim(fgl)[2]-1

glass <- levels(fgl$type)
glass
color.lst <- as.numeric(fgl$type)
stars(fgl[,1:p],full=FALSE,col.stars=color.lst)

levels.n <- length(glass)
means <- matrix(0,levels.n,p)
for (i in 1:levels.n){
	I <- fgl$type == glass[i]
	means[i,] <- colMeans(fgl[I,1:p])}
means <- data.frame(means,row.names =glass)
stars(means,full=FALSE,key.loc=c(3.2,3.2),key.labels =names(fgl)[1:p],byrow=T,col.stars=1:6)

Z <- scale(fgl[,1:p])

########### pair-wise comparisons: cycle x.var from 1 to 9 ###########
x.var <- 1
X <- rep(Z[,x.var],p)
Y <- as.vector(Z[,1:p])
G <- rep(fgl$type,p)
v.names <- names(fgl)[1:p]
variable <- rep(v.names,each=dim(Z)[1])
xyplot(Y~X|variable,groups=G,pch=16,xlab=names(fgl)[x.var])

R <- t(Z)%*%Z/(n-1)
R
###################### dimension reduction ###############

e.obj <- eigen(R)
e.ve <- e.obj$vectors
e.va <- e.obj$values

############## best approximations of R ################

R.approx <- e.ve[,1]%*%t(e.ve[,1])*e.va[1]
sum(abs(R))
sum(abs(R-R.approx))

for (v in 2:p){
	R.approx <- R.approx + e.ve[,v]%*%t(e.ve[,v])*e.va[v]
	print(sum(abs(R-R.approx)))}

#################### "variance explained by each eigenvector/axis #########
cbind(e.va/sum(e.va),cumsum(e.va),100*cumsum(e.va)/p)

################# apply dimension reduction ###############
P.matrix <- cbind(e.ve[,1],e.ve[,2],e.ve[,3])
P.matrix <- cbind(e.ve[,1]*e.va[1],e.ve[,2]*e.va[2],e.ve[,3]*e.va[3])

proj.1 <- Z%*%P.matrix
x <- proj.1[,1]
y <- proj.1[,2]

xyplot(y~x,groups=fgl$type,pch=16)

################## look at all 3 combinations

n <- dim(Z)[1]
x <- c(proj.1[,1],proj.1[,1],proj.1[,2])
y <- c(proj.1[,2],proj.1[,3],proj.1[,3])
G <- rep(fgl$type,3)
variables <- rep(c("1 vs 2","1 vs 3","2 vs 3"),each=n)
dim(cbind(x,y,G,variables))
xyplot(y~x|variables,groups=G,pch=16)