Why do the factor loadings from R's fa() and factanal() functions differ when using the promax rotation?

Question

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.

Answer 1

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