In a related question, I have asked why there are differences between stats::varimax and GPArotation::Varimax, both of which psych::principal calls, depending on the option set for rotate =.
The differences between these two (see other question) account for some, but not all of the differences from psych::principal.
It appears that these differences somehow get exacerbated by psych::principal.
(I have a simple theory why, and I'd like to get that confirmed).
library(GPArotation)
library(psych)
data("Thurstone")
principal.unrotated <- principal(r = Thurstone, nfactors = 4, rotate = "none") # find unrotated PCs first
loa <- unclass(principal.unrotated$loadings)
# just to compare that the rot.mat is correct
varimax.stats <- stats::varimax(x = loa, normalize = TRUE)
varimax.GPA <- GPArotation::Varimax(L = loa, normalize = TRUE)
# notice we're here NOT interested in the difference between stats and GPA, that's the other question
diff.from.rot.meth <- unclass(varimax.stats$loadings - varimax.GPA$loadings) # very small differences, see this question: https://stackoverflow.com/questions/32350891/why-are-there-differences-between-gparotationvarimax-and-statsvarimax
mean(abs(diff.from.rot.meth))
#> [1] 8.036863e-05
principal.varimax.stats <- principal(r = Thurstone, nfactors = 4, rotate = "varimax")
principal.Varimax.GPA <- principal(r = Thurstone, nfactors = 4, rotate = "Varimax")
diff.from.princ <- principal.Varimax.GPA$rot.mat - principal.varimax.stats$rot.mat # quite a substantial change, because Theta is NOT a rotmat, that makes sense
mean(abs(diff.from.princ))
#> [1] 0.021233
mean(abs(diff.from.rot.meth)) - mean(abs(diff.from.princ)) # principal has MUCH bigger differences
#> [1] -0.02115263
This seems way too big for just a floating point artefact or something.
My hypothesis would be that the (additional) difference stems from the fact that GPArotation::Varimax defaults to (Kaiser) normalize == FALSE, whereas **stats::varimax defaults to (Kaiser) normalize == TRUE, which cannot be set differently from within `principal::psych``.
stats::varimax manual:
## varimax with normalize = TRUE is the default
GPArotation::Varimax / GPArotation::GPForth manual:
The argument normalize gives an indication of if and how any normalization should be done before rotation, and then undone after rotation. If normalize is FALSE (the default) no normalization is done. If normalize is TRUE then Kaiser normalization is done. (So squared row entries of normalized A sum to 1.0. This is sometimes called Horst normalization.)
Also, they psych::Kaiser manual warns:
The GPArotation package does not (by default) normalize, nor does the fa function. Then, to make it more confusing, varimax in stats does,Varimax in GPArotation does not.
Can anyone confirm that the difference is, in fact, accounted for by the normalization options?
The problem of different loadings occurs because of different accuracy in the procedures (of course, since in psych::principal the normalize option is not evaluated, all the other procedures must used with that option switched to TRUE). While the accuracy in stats::varimax and GPArotation::Varimax can be configured (parameter eps) this is ignored in psych::principal, and seems to be implicitely equal to stats::varimax with eps=1e-5 .
If we increase accuracy in stats::varimax and GPArotations::Varimax to eps=1e-15 then we get the same result (at least up to 8 digits) as in a concurring implementation in my MatMate-program which I've tested to be very accurately equal to the SPSS computations too.
The missing handling of the explicite option eps in psych::principal seems to be a bug, and its poor implicite default value is surely unsatisfying.
Interestingly, the GPArotation::Varimax needs very many rotations with eps=1e-15, see the output at the end; so either there is another internal procedure implemented or the eps-parameter is evaluated differently in the decision when the iterations can be stopped. One example, see at the end of this answer, might suggest such an effect.
Below the protocol of the comparision; only the first two rows of the loadings are displayed.
The first three computations with low accuracy (eps=1e-5)
all procedures give equal or comparable results,
except GPArotation, which is already near the exact answer
--------------------------------
> principal(r = Thurstone, nfactors = 4, rotate = "varimax")$loadings *1000000
> stats::varimax(x = loa, normalize = TRUE,eps=1e-5)$loadings *1000000
> GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-5)$loadings *1000000
The second three computations with (attempted) high accuracy (eps=1e-15)
all procedures except psych::principal give equal results, and they
agree also with external crosscheck using MatMate
--------------------------------
> principal(r = Thurstone, nfactors = 4, rotate = "varimax",eps=1e-15)$loadings*1000000
> stats::varimax(x = loa, normalize = TRUE,eps=1e-15)$loadings*1000000
> GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-10)$loadings*1000000
The attempts to little/default accurate results with large eps or without
# ===== Numerical documentation (only first two rows are displayed)================
> principal(r = Thurstone, nfactors = 4, rotate = "varimax")$loadings*1000000
RC1 RC2 RC3 RC4
Sentences 871655.72 216638.46 198427.07 175202.57
Vocabulary 855609.28 294166.99 153181.45 180525.99
> stats::varimax(x = loa, normalize = TRUE,eps=1e-5)$loadings *1000000
PC1 PC2 PC3 PC4
Sentences 871655.72 216638.46 198427.07 175202.57
Vocabulary 855609.28 294166.99 153181.45 180525.99
> GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-5)$loadings *1000000
PC1 PC2 PC3 PC4
Sentences 871717.3 216618.7 198176.3 175204.47
Vocabulary 855663.1 294146.3 152930.7 180517.21
# =============================================================================
Now the attempts to more accurate results with smaller eps
> principal(r = Thurstone, nfactors = 4, rotate = "varimax",eps=1e-15)$loadings *1000000
RC1 RC2 RC3 RC4
Sentences 871655.72 216638.46 198427.07 175202.57
Vocabulary 855609.28 294166.99 153181.45 180525.99
> stats::varimax(x = loa, normalize = TRUE,eps=1e-15)$loadings *1000000
PC1 PC2 PC3 PC4
Sentences 871716.83 216619.69 198174.31 175207.86
Vocabulary 855662.58 294147.47 152928.77 180519.37
> GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-10)$loadings *1000000
PC1 PC2 PC3 PC4
Sentences 871716.8 216619.7 198174.31 175207.86
Vocabulary 855662.6 294147.5 152928.77 180519.37
Warnmeldung:
In GPForth(L, Tmat = Tmat, method = "varimax", normalize = normalize, :
convergence not obtained in GPForth. 1000 iterations used.
# Result by MatMate: --------------------------------------------------------
lad = cholesky(Thurstone)
pc = rot(lad,"pca")
pc4 = pc[*,1..4] // arrive at the first four pc's
t = gettrans( normzl(pc4),"varimax") // get rotation-matrix for row-normalized pc's
vmx = pc4 * t // rotate pc4 by rotation-matrix
display = vmx * 1000000
PC1 PC2 PC3 PC4
Sentences 871716.83 216619.68 198174.31 175207.87
Vocabulary 855662.58 294147.46 152928.77 180519.37
# ===============================================================================>
A much better matching of the results by stats::varimax and GPArotation::Varimax could be achieved by setting eps to 1e-12 in stat:: and to 1e-6 in GPArotation , where one value is the square of the other. We get then
> GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-6)$loadings*1000000
PC1 PC2 PC3 PC4
Sentences 871716.8 216619.8 198174.49 175207.63
Vocabulary 855662.5 294147.6 152928.94 180519.21
> stats::varimax(x = loa, normalize = TRUE,eps=1e-12)$loadings*1000000
PC1 PC2 PC3 PC4
Sentences 871716.80 216619.74 198174.40 175207.85
Vocabulary 855662.55 294147.52 152928.86 180519.36