## extra - BIPLOT


biplot(prcomp(X,scale=FALSE))
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")





color_biplot<-function(data,
                       axes=1:2,
                       cex=1,
                       ind.supp=-c(1:nrow(data)),
                       var.supp=-c(1:ncol(data)),
                       col.ind=rep("black",nrow(data)),
                       col.var=rep("red",ncol(data))){
  data<-scale(data,scale=FALSE)
  M<-nrow(data)
  q<-ncol(data)
  X<-data[-ind.supp,-var.supp]
  N <-nrow(X)
  p <-ncol(X)
  X.svd <- svd(X)
  X.pca<-prcomp(X)
  U <- X.svd$u
  D <- diag(X.svd$d)
  V<- X.svd$v
  scores <- U %*% D   ### = X.pca$x
  sdev<-diag(D)/sqrt(N-1) ### PC's standard deviations = X.pca$sdev
  V.aux<-V
  if (q>p){
    for (j in var.supp){
      V.aux<- rbind(V.aux,cor(data[,j],scores)/sdev*sd(data[,j]))
    }
  }
  G<-data %*% V.aux %*% solve(D)*sqrt(N-1)
  # G=U*sqrt(N-1)     
  H=V.aux %*% D  / sqrt(N-1)
  plot(G,cex=.001,asp=TRUE,xlab="PC1",ylab="PC2")
  abline(h=0,lty=2,col="gray")
  abline(v=0,lty=2,col="gray")
  text(G[,1],G[,2],col=col.ind,cex=cex)
  arrows(rep(0,p),rep(0,p),H[,1],H[,2],col=col.var,length=.1)
  text(H[,1]*1.3,H[,2]*1.3,pch=4,col=col.var,cex=1,labels=colnames(X))
}

### exemplo de utilizacao
color_biplot(iris[-5],col.ind = rep(c("red","blue","darkgreen"),c(50,50,50)))

