There are several implementations of factor analysis in Python. The following code returns the loadings obtained with three of these methods:
import seaborn as sns
from sklearn.preprocessing import StandardScaler
# THE DATA:
iris = sns.load_dataset('iris')
scaler = StandardScaler()
iris = scaler.fit_transform(iris.iloc[:,:-1])
NUM_COMPONENTS = 1
# METHOD 1:
from sklearn.decomposition import FactorAnalysis
fa1 = FactorAnalysis(n_components=NUM_COMPONENTS, svd_method='lapack', rotation=None)
_ = fa1.fit_transform(iris).flatten()
fa1.components_.T
# METHOD 2:
from factanal.wrapper import factanal
fa2 = factanal(iris, factors=NUM_COMPONENTS, scores='regression', rotation='none', verbose=True, return_dict=True)
fa2['loadings']
# METHOD 3:
from factor_analyzer.factor_analyzer import FactorAnalyzer
fa3 = FactorAnalyzer(n_factors=NUM_COMPONENTS, rotation=None, method='ml')
_ = fa3.fit_transform(iris)
fa3.loadings_
With one factor (i.e. NUM_COMPONENTS = 1), all three methods return the same loadings (ignoring sign differences). However, increasing the number of factors leads to different results. Method 2 throws an error because the degrees of freedom is negative, and the loadings for Methods 1 and 3 are very different.
Does anyone have any insight into why this is the case? Does it have to do with the negative degrees of freedom, or is it simply a difference in the optimisation procedure?
Note: If you replace the iris data with random data (np.random.normal(loc=0, scale=1, size=(150,5) using np.random.seed(12345)), Methods 2 and 3 return the same loadings when NUM_COMPONENTS is 1 or 2 (+ve d.o.f), but these are different to those returned by Method 1. So it seems like the sklearn implementation sometimes returns different results to the other two methods, even with one factor. Why is this?