Suppose I have two matrices of PCA loadings loa.orig, and loa.rot, and I know that loa.rot is a rotation (by-hand or otherwise) of loa.orig.
(loa.orig might also have been already orthogonally rotated by varimax or something, but I don't think that matters).
I know want to know the angles by which loa.orig has been rotated to arrive at loa.rot.
I understand from this comment on another question that "rotations don't commute", and therefore, the order of pair-wise (plane-wise) rotations also matters.
So to reproduce loa.rot from loa.orig I would need to know a series of necessary rotations, ideally in the order given by rots in the below.
Here's an MWE:
library(psych) # this allows manual rotation
# suppose I have some ORIGINAL loadings matrix, from a principal components analysis, with three retained components
loa.orig <- cbind(c(0.6101496, 0.7114088, 0.3356003, 0.7318809, 0.5980133, 0.4102817, 0.7059148, 0.6080662, 0.5089014, 0.587025, 0.6166816, 0.6728603, 0.7482675, 0.5409658, 0.6415472, 0.3655053, 0.6313868), c(-0.205317, 0.3273207, 0.7551585, -0.1981179, -0.423377, -0.07281187, -0.04180098, 0.5003459, -0.504371, 0.1942334, -0.3285095, 0.5221494, 0.1850734, -0.2993066, -0.08715662, -0.02191772, -0.2002428), c(-0.4692407, 0.1581682, -0.04574932, -0.1189175, 0.2449018, -0.5283772, 0.02826476, 0.1703277, 0.2305158, 0.2135566, -0.2783354, -0.05187637, -0.104919, 0.5054129, -0.2403471, 0.5380329, -0.07999642))
# I then rotate 1-2 by 90°, and 1-3 by 45°
loa.rot <- factor.rotate(f = loa.orig, angle = 90, col1 = 1, col2 = 2)
loa.rot <- factor.rotate(f = loa.rot, angle = 45, col1 = 1, col2 = 3)
# predictably, loa.rot and loa.orig are now different
any(loa.rot == loa.orig) # are any of them the same?
Obviously, in this case, I know the angles and order, but let's assume I don't. Also, let's assume that in the real use case there may be many components retained and rotated, not just three.
I am a bit unsure what would be a conventional way to report the order of component-pair-wise (plane-wise) rotation angles, but I'd imagine the list of possible combinations (~~not permutations~~) should do.
combs <- combn(x = ncol(loa.orig), m = 2, simplify = TRUE) # find all possible combinations of factors
rots <- data.frame(t(combs), stringsAsFactors = FALSE) # transpose
rots # these rows give the *order* in which the rotations should be done
rots gives these permutations.
Would be great to know how to arrive at loa.rot from loa.orig, rotating component pairs as given by the rows in rots.
UPDATE: Attempt based on below answer
Based on the below answer, I've tried to put together a function and test it with a varimax rotation and a real dataset.
(No reason in particular for varimax – I just wanted some rotation where we actually do not know the angles.).
I then test whether I can actually re-create the varimax rotation from the vanilla loadings, using the extracted angles.
library(psych)
data("Harman74.cor") # notice the correlation matrix is called "cov", but doc says its a cor matrix
vanilla <- principal(r = Harman74.cor$cov, nfactors = 4, rotate = "none", )$loadings # this is unrotated
class(vanilla) <- NULL # print methods only causes confusion
varimax <- principal(r = Harman74.cor$cov, nfactors = 4, rotate = "varimax")$loadings # this is rotated
class(varimax) <- NULL # print methods only causes confusion
find.rot.instr <- function(original, rotated) {
# original <- vanilla$loadings # testing
# rotated <- varimax$loadings # testing
getAngle <- function(A, B) acos(sum(A*B) / (norm(A, "F") * norm(B, "F"))) * 180/pi
rots <- combn(x = ncol(original), m = 2, simplify = FALSE) # find all possible combinations of factor pairs
tmp <- original
angles <- sapply(rots, function(cols) {
angle <- getAngle(tmp[, cols], rotated[, cols])
tmp <<- factor.rotate(tmp, angle = angle, col1 = cols[1], col2 = cols[2])
return(angle)
})
return(angles)
}
vanilla.to.varimax.instr <- find.rot.instr(original = vanilla, rotated = varimax) # these are the angles we would need to transform in this order
rots <- combn(x = ncol(vanilla), m = 2, simplify = FALSE) # find all possible combinations of factor pairs
# this is again, because above is in function
# now let's implement the extracted "recipe"
varimax.recreated <- vanilla # start with original loadings
varimax.recreated == vanilla # confirm that it IS the same
for (i in 1:length(rots)) { # loop over all combinations, starting from the top
varimax.recreated <- factor.rotate(f = varimax.recreated, angle = vanilla.to.varimax.instr[i], col1 = rots[[i]][1], col2 = rots[[i]][2])
}
varimax == varimax.recreated # test whether they are the same
varimax - varimax.recreated # are the close?
Unfortunately, they are not the same, not even similar :(
> varimax == varimax.recreated # test whether they are the same
PC1 PC3 PC2 PC4
VisualPerception FALSE FALSE FALSE FALSE
Cubes FALSE FALSE FALSE FALSE
PaperFormBoard FALSE FALSE FALSE FALSE
Flags FALSE FALSE FALSE FALSE
GeneralInformation FALSE FALSE FALSE FALSE
PargraphComprehension FALSE FALSE FALSE FALSE
SentenceCompletion FALSE FALSE FALSE FALSE
WordClassification FALSE FALSE FALSE FALSE
WordMeaning FALSE FALSE FALSE FALSE
Addition FALSE FALSE FALSE FALSE
Code FALSE FALSE FALSE FALSE
CountingDots FALSE FALSE FALSE FALSE
StraightCurvedCapitals FALSE FALSE FALSE FALSE
WordRecognition FALSE FALSE FALSE FALSE
NumberRecognition FALSE FALSE FALSE FALSE
FigureRecognition FALSE FALSE FALSE FALSE
ObjectNumber FALSE FALSE FALSE FALSE
NumberFigure FALSE FALSE FALSE FALSE
FigureWord FALSE FALSE FALSE FALSE
Deduction FALSE FALSE FALSE FALSE
NumericalPuzzles FALSE FALSE FALSE FALSE
ProblemReasoning FALSE FALSE FALSE FALSE
SeriesCompletion FALSE FALSE FALSE FALSE
ArithmeticProblems FALSE FALSE FALSE FALSE
> varimax - varimax.recreated # are the close?
PC1 PC3 PC2 PC4
VisualPerception 0.2975463 1.06789735 0.467850675 0.7740766
Cubes 0.2317711 0.91086618 0.361004861 0.4366521
PaperFormBoard 0.1840995 0.98694002 0.369663215 0.5496151
Flags 0.4158185 0.82820078 0.439876777 0.5312143
GeneralInformation 0.8807097 -0.33385999 0.428455899 0.7537385
PargraphComprehension 0.7604679 -0.30162120 0.389727192 0.8329341
SentenceCompletion 0.9682664 -0.39302764 0.445263121 0.6673116
WordClassification 0.7714312 0.03747430 0.460461099 0.7643221
WordMeaning 0.8010876 -0.35125832 0.396077591 0.8201986
Addition 0.4236932 -0.32573100 0.204307400 0.6380764
Code 0.1654224 -0.01757153 0.194533996 0.9777764
CountingDots 0.3585004 0.28032822 0.301148474 0.5929926
StraightCurvedCapitals 0.5313385 0.55251701 0.452293566 0.6859854
WordRecognition -0.3157408 -0.13019630 -0.034647588 1.1235253
NumberRecognition -0.4221889 0.10729098 -0.035324356 1.0963785
FigureRecognition -0.3213392 0.76012989 0.158748259 1.1327322
ObjectNumber -0.3234966 -0.02363732 -0.007830001 1.1804147
NumberFigure -0.2033601 0.59238705 0.170467459 1.0831672
FigureWord -0.0788080 0.35303097 0.154132395 0.9097971
Deduction 0.3423495 0.41210812 0.363022937 0.9181519
NumericalPuzzles 0.3573858 0.57718626 0.393958036 0.8206304
ProblemReasoning 0.3430690 0.39082641 0.358095577 0.9133117
SeriesCompletion 0.4933886 0.56821932 0.465602192 0.9062039
ArithmeticProblems 0.4835965 -0.03474482 0.332889805 0.9364874
So clearly, I'm making a mistake.
I think the problem is easier if you do the matrix algebra (I am doing this for orthogonal rotations, not oblique transformations, but the logic is the same.)
First, notice that any rotation, t, produces a rotated components matrix c such that c = Ct where C is the unrotated PCA solution.
Then, since C'C is the original correlation matrix, R, we can solve for t by premultiplying both sides by C' and then by the inverse of R. This leads to
t = C' R^-1 c
For your example, let
R <- Harman74.cor$cov
P <- principal(Harman74.cor$cov,nfactors=4,rotate="none")
p <- principal(Harman74.cor$cov,nfactors=4) #the default is to do varimax
rotation <- t(P$loadings) %*% solve(R) %*% p$loadings
Then, to see if this correct
P$loadings %*% rotation - p$loadings
In addition, the coming release of psych now reports the rotation matrix, so we can compare rotation to p$rot.mat
round(p$rot.mat,2)
round(rotation,2)
These differ in the order of the components, but that is because psych reorders the components after rotation to reflect the size of the rotated component sum of squares.