Last active
June 9, 2016 16:06
-
-
Save rmaia/be85ff1f0a89c65207870db3e1308056 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment