eigensinplot.R

# Plotting eigenvectors from PCA onto scatterplots 

sigmat <- matrix(c(1,-.9,-.9,1), nrow=2)
cent <- c(8,3)

dat <- MASS::mvrnorm(n=100, mu=cent, Sigma=sigmat)

eig <- eigen(cov(dat))

plot(dat, xlim=cent[1]+c(-3,3)*sigmat[1,1], 
    ylim=cent[2]+c(-3,3)*sigmat[2,2], pch=19, col=grey(0.8))



# important concepts to remember:

# Eigenvectors: 
# - Angles (in radians) between PC and original variables
# - cells representing loadings of variable i on PCj 

# Eigenvalues: 
# - Are variances of the new variables
# - var = SD^2
# - If data are multivariate normal, it'll have a range of approximately 4*SD
# - Therefore ~95% of the data should be between SD*c(-2,2) of the mean


scaleig <- 2*sqrt(eig$value)
meaneig <- colMeans(dat)

# For PCA1, we'll use the first eigenvalue and the first column of eigenvectors
# add mean to center on data

segments(x0= meaneig[1] - scaleig[1]*eig$vector[1,1],
         y0= meaneig[2] - scaleig[1]*eig$vector[2,1],
         x1= meaneig[1] + scaleig[1]*eig$vector[1,1],
         y1= meaneig[2] + scaleig[1]*eig$vector[2,1],
         lwd=3)
         
# For PCA2, we'll use the second eigenvalue and the second column of eigenvectors
# add mean to center on data

segments(x0= meaneig[1] - scaleig[2]*eig$vector[1,2],
         y0= meaneig[2] - scaleig[2]*eig$vector[2,2],
         x1= meaneig[1] + scaleig[2]*eig$vector[1,2],
         y1= meaneig[2] + scaleig[2]*eig$vector[2,2],
         lwd=3)

# validating
# calculated different ways (got from the internet)

eigScl  <- eig$vector %*% diag(sqrt(eig$value))  # scale eigenvectors to length = square-root
xMat    <- rbind(meaneig[1] + eigScl[1, ], meaneig[1] - eigScl[1, ])
yMat    <- rbind(meaneig[2] + eigScl[2, ], meaneig[2] - eigScl[2, ])
matlines(xMat, yMat, lty=1, lwd=2.5, col="#ec7014")

car::ellipse(center=meaneig, shape=sigmat, radius=2, col="#238443", lty=2, center.cex=0)