After doing PCA, the first component describes the largest part of variability. This is important e.g. in study of body measurements where it is commonly known (Jolliffe, 2002) that PC1 axis captures size variation. My question is whether PCA scores after varimax rotation retain the same properties or are they different as mentioned in this topic?

Since I need PCA scores for further statistical analyses I am wondering if varimax is needed and does it in fact disrupt the representation of real sample variability so that individual scores on rotated axes are uninformative or lead to miss-interpretation of reality?

Also could someone suggest some other references on this topic?

Workflows in R:

1. PCA (FactoMineR or prcomp) -> Extract individual scores -> Enter scores in the lm
2. PCA (FactoMiner or prcomp) -> Varimax on loadings matrix -> calculate the individual scores -> enter scores in the lm
3. FA (psych, varimax and pca extraction method) -> extract individual scores -> Enter scores in the lm

Now, without rotation (1.) percentages of explained variability are i.e. 29.32, 5.6, 3.2, on the first three axes. 2. and 3. solutions yield similar percentages on the first three factors i.e. 12.2, 12.1, 8.2. Off course 1. solution tends to push all high variable loadings on the first axis, while 2. and 3. tend to distribute loadings between axes (which is the reason for rotation). I wanted to know if these three workflows are essential the same since individual scores are different on rotated vs. unrotated axes?