# r code adapted from Antoni Parellada╒s post on: http://stats.stackexchange.com/questions/2691/making-sense-of-principal-component-analysis-eigenvectors-eigenvalues

mydata<-read.table("mydata.csv", header=T, sep=",", row.names=1) #loads the .csv file
d=data.frame(x=mydata[,1], y=mydata[,2])
mean_dev_data <- as.data.frame(d[,1:2])
mean_dev_data
X          <- as.matrix(mean_dev_data) 
S             <- cov(X)                               # Covariance matrix
lambda        <- eigen(S)$values                      # Eigenvalues
lambda_matrix <- diag(2)*eigen(S)$values              # Eigenvalues matrix
e_vectors     <- eigen(S)$vectors                     # Eigenvectors
colnames(e_vectors) <- c("PC1","PC2")
e_vectors[,1] <- e_vectors[,1]*-1
score_matrix <-  X %*% e_vectors          # Identical to the often found operation:
t(t(e_vectors) %*% t(X))
# For scores:
  prcomp(d, center=T, scale=T) # or prcomp(mean_dev_data)
    # or...
    princomp(mean_dev_data)$scores

# and for eigenvectors:
  prcomp(mean_dev_data)$rotation
    # or...
    princomp(mean_dev_data)$loadings

# and for eigenvalues:
  prcomp(mean_dev_data)$sdev^2
    # or...
    princomp(covmat=C)$sd^2

plot(x~y, main="Scaled Scores on 2 Docuscope LATs", pch=19, data=mean_dev_data, xlab="Abstract Concepts", ylab="First Person", xlim=c(-.8,.8), ylim=c(-.8,.8))
abline(v=0, h=0, col="light gray")


plot(x~y, pch=19, main="Scaled Scores with Eigenvector Directions 1", data=mean_dev_data, xlab="Abstract Concepts", ylab="First Person", xlim=c(-.8,.8), ylim=c(-.8,.8))
abline(v=0, h=0, col="light gray")
#with(mean_dev_data, text(x~y, pos=3)) #item labels
abline(0, e_vectors[2,1]/e_vectors[1,1], col="red")
#arrows(x0=0, x1=lambda[1]*e_vectors[1,1],0, y1=lambda[1]*e_vectors[2,1], col="red",  lwd=2)
#arrows(x0=0, x1=lambda[2]*e_vectors[1,2], y0=0, y1=lambda[2]*e_vectors[2,2], col="black", lwd=2)
line1 <- c(0, e_vectors[2,1]/e_vectors[1,1])
perp.segment.coord <- function(x0, y0, line1){
  a <- line1[1]  # intercept
  b <- line1[2]  # slope

  x1 <- (x0 + b*y0 - a*b)/(1 + b^2)
  y1 <- a + b*x1
  list(x0=x0, y0=y0, x1=x1, y1=y1)
}
ss <- perp.segment.coord(mean_dev_data$y, mean_dev_data$x, line1)
segments(x0=ss$x0, x1=ss$x1, y0=ss$y0, y1=ss$y1, col="dark gray")

plot(x~y, pch=19, main="Scaled Scores with Eigenvector Directions 2", data=mean_dev_data, xlab="Abstract Concepts", ylab="First Person", xlim=c(-.8,.8), ylim=c(-.8,.8))
abline(v=0, h=0, col="dark gray")
#with(mean_dev_data, text(x~y, pos=3)) #item labels
abline(0, e_vectors[2,2]/e_vectors[1,2], col="black")
line1a <- c(0, e_vectors[2,2]/e_vectors[1,2]) 
ssa    <- perp.segment.coord(mean_dev_data$y, mean_dev_data$x, line1a) 
segments(x0=ssa$x0, x1=ssa$x1, y0=ssa$y0, y1=ssa$y1, col="light gray")

score_matrix[,1] = score_matrix[,1]*-1 #flip signs of x axis
score_matrix[,2]=score_matrix[,2]*-1 # flip signs of y axis
plot (score_matrix, pch=19, main="Data Replotted on New Bases (PC1 and PC2)", xlim=c(-.8,.8), ylim=c(-.8,.8))
abline(a=0, b=0, col="red")

line1 <- c(0, 0)
perp.segment.coord <- function(x0, y0, line1){
  a <- line1[1]  # intercept
  b <- line1[2]  # slope

  x1 <- (x0 + b*y0 - a*b)/(1 + b^2)
  y1 <- a + b*x1
  list(x0=x0, y0=y0, x1=x1, y1=y1)
}
ss <- perp.segment.coord(score_matrix[,1], score_matrix[,2], line1)
segments(x0=ss$x0, x1=ss$x1, y0=ss$y0, y1=ss$y1, col="dark gray")

score_matrix0<-cbind(score_matrix[,1],0)
plot (score_matrix0, pch=19, main="Reduction to PC1", xlim=c(-.8,.8), ylim=c(-.8,.8))
abline(a=0, b=0, col="red")