# rotate "perfekt" data by -60 degrees library("mvtnorm") z <- rmvnorm(1000,sigma=matrix(c(5,0,0,1),nrow=2)) RotMat <- matrix(c(cos(pi/3),-sin(pi/3),sin(pi/3),cos(pi/3)),nrow=2) x <- z %*% RotMat plot(z,xlim=c(-7,7),ylim=c(-7,7)) points(x,col="red") abline(b=tan(-pi/3),a=0) pca <- prcomp(x) points(pca$x,col="yellow") names(pca) pca ( pca$sdev )^2 RotMat %*% pca$rotation t( pca$rotation ) %*% pca$rotation cov(z) t( pca$rotation ) %*% cov(x) %*% pca$rotation ####################################################### # Height shoe weight data: shsw <-read.table("HeightShoeWeight.txt",header=TRUE) attach(shsw) head(shsw) hsw <- shsw[,2:4] head(hsw) hsw.pca <- prcomp(hsw,scale=TRUE) hsw.pca # proportion of variance explained by first component hsw.pca$sdev[1]^2/sum(hsw.pca$sdev^2) fm.col <- character() fm.col[sex==0] <- "blue" fm.col[sex==1] <- "red" sqrt( length(sex)-1 ) # = 15 plot(hsw.pca$x,ylim=c(-3,6)) plot(hsw.pca$x,ylim=c(-3,6),col=fm.col) biplot(hsw.pca,scale=0) biplot(hsw.pca,scale=1) biplot(hsw.pca,scale=1,xlabs=sex) # Relative contribution of the different variables # to the first principal component: v <- t(hsw.pca$rotation[,1]) abs(v)/sum(abs(v)) ####################################################### # Countries data countries <- read.table("Countries.txt",header=TRUE) cntr.pca <- prcomp(countries,scale=TRUE); cntr.pca plot(cntr.pca$x) biplot(cntr.pca,scale=0) biplot(cntr.pca,scale=1) ####################################################### # Height shoe weight data hsw.pca$sdev^2/sum( (hsw.pca$sdev)^2 ) cumsum( hsw.pca$sdev^2/sum( (hsw.pca$sdev)^2 ) ) screeplot( hsw.pca, type="lines") hsw.pca$sdev^2/sum( (hsw.pca$sdev)^2 ) p<-length(hsw.pca$sdev) L<-matrix(ncol=p) for (i in 1:p) { L[i] <- round(1/p*sum(1/seq(from=i, to=p)),2)} L ####################################################### # EWU data ewu <- read.table("EWU.txt",header=TRUE) ewu1 <- ewu[,2:5] ewu.pca <- prcomp(ewu1, scale=TRUE) biplot(ewu.pca,scale=0,xlabs=ewu$Staat) biplot(ewu.pca,scale=1,xlabs=ewu$Staat)