I tried to use factanal() from base R and fa() from the psych package to perform a factor analysis on data from a questionnaire with same response scale for each question.
Why do I obtain different factor loadings from the two functions even with having Maximum likelihood set for both methods when using the "promax" rotation? This is a reproducible example:
# Load the psych package for
# data analysis and visualization
library(psych)
library(GPArotation)
# Load the mtcars dataset
data(mtcars)
# Perform factor analysis on the mtcars dataset
factanal <- factanal(mtcars, factors=3, rotation="promax")
fa <- fa(mtcars, nfactors = 3, fm="ml", rotate="promax")
# Print the results
head(round(factanal$loadings, 2))
Factor1 Factor2 Factor3
mpg 0.54 0.16 -0.43
cyl -0.53 -0.58 0.08
disp -0.65 -0.32 0.19
hp -0.11 -0.48 0.47
drat 0.83 0.12 0.10
wt -0.71 0.18 0.54
head(round(fa$loadings, 2))
ML2 ML1 ML3
mpg 0.52 0.17 -0.43
cyl -0.51 -0.59 0.08
disp -0.63 -0.33 0.19
hp -0.09 -0.49 0.48
drat 0.83 0.13 0.10
wt -0.70 0.16 0.54
In R documentation, they recommend to use "oblique.Scores=TRUE", but that does not change the factor loadings:
fa2 <- fa(mtcars, nfactors = 3, fm="ml", rotate="promax", oblique.scores=T, scores="tenBerge")
head(round(fa2$loadings, 2))
ML2 ML1 ML3
mpg 0.52 0.17 -0.43
cyl -0.51 -0.59 0.08
disp -0.63 -0.33 0.19
hp -0.09 -0.49 0.48
drat 0.83 0.13 0.10
wt -0.70 0.16 0.54
While the question has already been asked here https://stackoverflow.com/questions/50573764/differences-between-fa-and-factanal-functions-in-r#:~:text=factanal%20performs%20a%20maximum%2Dlikelihood,least%20square%20regressions%20(OLR). the answer is not completely true, as in the article they link to they say "same results" but it is same results only when we use "varimax" rotation.
You can get the same rotation as stats::promax, used by factanal, by using Promax (note the upper case P). The help page of ?Promax explains:
Promax is a very direct adaptation of the stats::promax function. The addition is that it will return the interfactor correlations as well as the loadings and rotation matrix.
I have not looked in the docs to see what promax, from the psych package, is doing (I'm sure it will be in there as they are very good) but a quick look at the source code for fa reveals that it returns a rotation of the normalised loadings:
Promax = {pro <- Promax(loadings,...) #Promax without Kaiser normalization
loadings <- pro$loadings
Phi <- pro$Phi
rot.mat <- pro$rotmat},
promax = {#pro <- stats::promax(loadings,...) #from stats
pro <- kaiser(loadings,rotate="Promax",...) #calling promax will now do the Kaiser normalization before doing Promax rotation
Code showing the equivalence:
library(psych)
data(mtcars)
m1 <- factanal(mtcars, factors=3, rotation="promax")
m2 <- fa(mtcars, nfactors = 3, fm="ml", rotate="Promax")
loadings(m2)[] - loadings(m1)[]
# ML2 ML1 ML3
# mpg -1.705875e-06 1.473799e-06 -3.059439e-07
# cyl 1.594591e-07 -8.884557e-07 -1.818443e-07
# disp 8.168498e-07 4.712348e-07 2.310251e-06
# hp 6.522325e-07 -1.931424e-06 -1.722668e-06
# drat 8.599102e-07 1.509280e-06 8.609054e-07
# wt 2.121750e-06 8.515195e-07 3.925057e-06
# qsec 1.092731e-06 1.119313e-06 2.633606e-06
# vs 4.688156e-07 -8.692953e-07 -6.737374e-07
# am -1.884347e-06 1.230842e-06 -1.514240e-06
# gear -1.691097e-06 -2.111156e-06 -5.330317e-06
# carb 1.288899e-06 -2.468629e-06 -2.641025e-06