I ran PCA on 25 variables and selected the top 7 PCs using
prcomp
.prc < prcomp(pollutions, center=T, scale=T, retx=T)
I have then done varimax rotation on those components.
varimax7 < varimax(prc$rotation[,1:7])
And now I wish to varimax rotate the PCArotated data (as it is not part of the varimax object – only the loadings matrix and the rotation matrix). I read that to do this you multiply the transpose of the rotation matrix by the transpose of the data so I would have done this:
newData < t(varimax7$rotmat) %*% t(prc$x[,1:7])
But that doesn’t make sense as the dimensions of the matrix transposes above are 7×7 and 7×16933 respectively and so I will be left with a matrix of only 7 rows, rather than 16933 rows… does anyone know what I am doing wrong here or what my final line should be? Do I just need to transpose back afterwards?
Answer
“Rotations” is an approach developed in factor analysis; there rotations (such as e.g. varimax) are applied to loadings, not to eigenvectors of the covariance matrix. Loadings are eigenvectors scaled by the square roots of the respective eigenvalues. After the varimax rotation, the loading vectors are not orthogonal anymore (even though the rotation is called “orthogonal”), so one cannot simply compute orthogonal projections of the data onto the rotated loading directions.
@FTusell’s answer assumes that varimax rotation is applied to the eigenvectors (not to loadings). This would be pretty unconventional. Please see my detailed account of PCA+varimax for details: Is PCA followed by a rotation (such as varimax) still PCA? Briefly, if we look at the SVD of the data matrix X=USV⊤, then to rotate the loadings means inserting RR⊤ for some rotation matrix R as follows: X=(UR)(R⊤SV⊤).
If rotation is applied to loadings (as it usually is), then there are at least three easy ways to compute varimaxrotated PCs in R :

They are readily available via function
psych::principal
(demonstrating that this is indeed the standard approach). Note that it returns standardized scores, i.e. all PCs have unit variance. 
One can manually use
varimax
function to rotate the loadings, and then use the new rotated loadings to obtain the scores; one needs to multiple the data with the transposed pseudoinverse of the rotated loadings (see formulas in this answer by @ttnphns). This will also yield standardized scores. 
One can use
varimax
function to rotate the loadings, and then use the$rotmat
rotation matrix to rotate the standardized scores obtained withprcomp
.
All three methods yield the same result:
irisX < iris[,1:4] # Iris data
ncomp < 2
pca_iris_rotated < psych::principal(irisX, rotate="varimax", nfactors=ncomp, scores=TRUE)
print(pca_iris_rotated$scores[1:5,]) # Scores returned by principal()
pca_iris < prcomp(irisX, center=T, scale=T)
rawLoadings < pca_iris$rotation[,1:ncomp] %*% diag(pca_iris$sdev, ncomp, ncomp)
rotatedLoadings < varimax(rawLoadings)$loadings
invLoadings < t(pracma::pinv(rotatedLoadings))
scores < scale(irisX) %*% invLoadings
print(scores[1:5,]) # Scores computed via rotated loadings
scores < scale(pca_iris$x[,1:2]) %*% varimax(rawLoadings)$rotmat
print(scores[1:5,]) # Scores computed via rotating the scores
This yields three identical outputs:
1 1.083475 0.9067262
2 1.377536 0.2648876
3 1.419832 0.1165198
4 1.471607 0.1474634
5 1.095296 1.0949536
Note: The varimax
function in R uses normalize = TRUE, eps = 1e5
parameters by default (see documentation). One might want to change these parameters (decrease the eps
tolerance and take care of Kaiser normalization) when comparing the results to other software such as SPSS. I thank @GottfriedHelms for bringing this to my attention. [Note: these parameters work when passed to the varimax
function, but do not work when passed to the psych::principal
function. This appears to be a bug that will be fixed.]
Attribution
Source : Link , Question Author : Scott , Answer Author : amoeba