Test of the hypothesis that the selected number of factors is sufficient

.pre-commit-config.yaml

@@ -3,7 +3,7 @@
 
 repos:
 -   repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.1.0
+    rev: v4.4.0
     hooks:
     -   id: trailing-whitespace
     -   id: end-of-file-fixer
@@ -14,27 +14,27 @@ repos:
     -   id: check-json
     -   id: debug-statements
     -   id: check-merge-conflict
--   repo: https://gitlab.com/pycqa/flake8
-    rev: 4.0.1
+-   repo: https://github.com/pycqa/flake8
+    rev: 6.0.0
     hooks:
     - id: flake8
       args: [--max-line-length=101]
       exclude: ^(factor_analyzer/__init__.py)
 -   repo: https://github.com/PyCQA/isort
-    rev: 5.10.1
+    rev: 5.12.0
     hooks:
     -   id: isort
         args: ["--profile", "black", "--filter-files"]
 -   repo: https://github.com/ikamensh/flynt/
-    rev: '0.76'
+    rev: '0.78'
     hooks:
     -   id: flynt
 -   repo: https://github.com/psf/black
-    rev: 22.3.0
+    rev: 23.3.0
     hooks:
     -   id: black
 -   repo: https://github.com/pycqa/pydocstyle
-    rev: 6.1.1
+    rev: 6.3.0
     hooks:
     -   id: pydocstyle
         args:

factor_analyzer/factor_analyzer.py

@@ -8,6 +8,7 @@
 """
 
 import warnings
+from typing import Tuple
 
 import numpy as np
 import pandas as pd
@@ -969,3 +970,36 @@ def get_factor_variance(self):
         check_is_fitted(self, "loadings_")
         loadings = self.loadings_.copy()
         return self._get_factor_variance(loadings)
+
+    def sufficiency(self, num_observations: int) -> Tuple[float, int, float]:
+        """
+        Perform the sufficiency test.
+
+        The test calculates statistics under the null hypothesis that
+        the selected number of factors is sufficient.
+
+        Parameters
+        ----------
+        num_observations: int
+            The number of observations in the input data that this factor
+            analyzer was fit using.
+
+        Returns
+        -------
+        statistic: float
+            The test statistic
+        degrees: int
+            The degrees of freedom
+        pvalue: float
+            The p-value of the test
+        """
+        nvar = self.corr_.shape[0]
+        degrees = ((nvar - self.n_factors) ** 2 - nvar - self.n_factors) // 2
+        obj = self._fit_ml_objective(
+            self.get_uniquenesses(), self.corr_, self.n_factors
+        )
+        statistic = (
+            num_observations - 1 - (2 * nvar + 5) / 6 - (2 * self.n_factors) / 3
+        ) * obj
+        pvalue = chi2.sf(statistic, df=degrees)
+        return statistic, degrees, pvalue

tests/expected/test01/sufficiency_ml_none_15_test01.csv

@@ -0,0 +1,16 @@
+n_factors,statistic,df,pvalue
+1,6158.69420371228,740,0
+2,3244.72903681528,701,4.44659081257122e-322
+3,1475.87556412021,663,8.80428369753958e-64
+4,1206.41930867432,626,3.36923786804043e-39
+5,960.29751580042,590,3.64376268923982e-20
+6,820.739852110472,555,1.36586692869634e-12
+7,720.71679931586,521,1.37467963009699e-08
+8,631.742866356578,488,1.14373979137026e-05
+9,558.423820827789,456,0.000715862260222117
+10,500.223805331245,425,0.00686044615605591
+11,446.535371500571,395,0.0373055704171928
+12,397.200698549131,366,0.125803730285426
+13,343.188514729509,338,0.411317090149109
+14,304.335701489558,311,0.595750888646432
+15,269.814900486032,285,0.732252710489368

tests/sufficiency_test.r

@@ -0,0 +1,15 @@
+data = read.csv("tests/data/test01.csv")
+
+factor_range = 1:15
+statistic = rep(0, length(factor_range))
+dof = rep(0, length(factor_range))
+pvalue = rep(0, length(factor_range))
+for (i in 1:length(factor_range)) {
+  res = factanal(data, factors=factor_range[i], rotation="none")
+  statistic[i] = res$STATISTIC
+  dof[i] = res$dof
+  pvalue[i] = res$PVAL
+}
+
+f = data.frame(n_factors=factor_range, statistic=statistic, df=dof, pvalue=pvalue)
+write.csv(f, "tests/expected/test01/sufficiency_ml_none_15_test01.csv", row.names=FALSE, quote=FALSE)

tests/test_factor_analyzer.py

@@ -24,7 +24,6 @@
 
 
 def test_calculate_bartlett_sphericity():  # noqa: D103
-
     path = "tests/data/test01.csv"
     data = pd.read_csv(path)
     s, p = calculate_bartlett_sphericity(data.values)
@@ -34,7 +33,6 @@ def test_calculate_bartlett_sphericity():  # noqa: D103
 
 
 def test_calculate_kmo():  # noqa: D103
-
     path = "tests/data/test02.csv"
     data = pd.read_csv(path)
 
@@ -83,7 +81,6 @@ def test_gridsearch():  # noqa: D103
 
 class TestFactorAnalyzer:  # noqa: D101
     def test_analyze_weights(self):  # noqa: D102
-
         data = pd.DataFrame(
             {
                 "A": [2, 4, 5, 6, 8, 9],
@@ -107,7 +104,6 @@ def test_analyze_weights(self):  # noqa: D102
         assert_array_almost_equal(expected_weights, fa.weights_)
 
     def test_analyze_impute_mean(self):  # noqa: D102
-
         data = pd.DataFrame(
             {
                 "A": [2, 4, 5, 6, 8, 9],
@@ -126,7 +122,6 @@ def test_analyze_impute_mean(self):  # noqa: D102
         assert_array_almost_equal(fa.corr_, expected_corr)
 
     def test_analyze_impute_median(self):  # noqa: D102
-
         data = pd.DataFrame(
             {
                 "A": [2, 4, 5, 6, 8, 9],
@@ -148,7 +143,6 @@ def test_analyze_impute_median(self):  # noqa: D102
         assert_array_almost_equal(expected_corr, fa.corr_)
 
     def test_analyze_impute_drop(self):  # noqa: D102
-
         data = pd.DataFrame(
             {
                 "A": [2, 4, 5, 6, 8, 9],
@@ -183,7 +177,6 @@ def test_analyze_rotation_value_error(self):  # noqa: D102
 
     @raises(ValueError)
     def test_analyze_infinite(self):  # noqa: D102
-
         data = pd.DataFrame(
             {
                 "A": [1.0, 0.4, 0.5],
@@ -215,7 +208,6 @@ def test_smc_is_r_squared(self):  # noqa: D102
         assert_array_almost_equal(smc_result, expected_r2)
 
     def test_factor_variance(self):  # noqa: D102
-
         path = "tests/data/test01.csv"
         data = pd.read_csv(path)
 
@@ -234,3 +226,25 @@ def test_factor_variance(self):  # noqa: D102
         proportional_variance = fa.get_factor_variance()[1]
 
         assert_array_almost_equal(proportional_variance_expected, proportional_variance)
+
+    def test_sufficiency(self):
+        path = "tests/data/test01.csv"
+        data = pd.read_csv(path)
+
+        # compute the sufficiency values for a number of factors
+        computed_values = []
+        for n_factors in range(1, 16):
+            fa = FactorAnalyzer(n_factors=n_factors, rotation=None, method="ml")
+            fa.fit(data)
+            computed_values.append([n_factors, *fa.sufficiency(data.shape[0])])
+
+        # create a data frame from the computed values
+        df_computed = pd.DataFrame(
+            computed_values, columns=["n_factors", "statistic", "df", "pvalue"]
+        )
+
+        # check against the values we expect
+        df_expected = pd.read_csv(
+            "tests/expected/test01/sufficiency_ml_none_15_test01.csv"
+        )
+        pd.testing.assert_frame_equal(df_computed, df_expected)

tests/test_utils.py

@@ -10,8 +10,7 @@
 import numpy as np
 import pandas as pd
 from nose.tools import eq_, raises
-from numpy.testing import assert_array_equal
-from pandas.util.testing import assert_almost_equal
+from numpy.testing import assert_array_equal, assert_almost_equal
 
 from factor_analyzer.utils import (
     commutation_matrix,
@@ -193,9 +192,9 @@ def test_partial_correlations():  # noqa: D103
     data = pd.DataFrame([[12, 14, 15], [24, 12, 52], [35, 12, 41], [23, 12, 42]])
 
     expected = [
-        [1.0, -0.730955, -0.50616],
-        [-0.730955, 1.0, -0.928701],
-        [-0.50616, -0.928701, 1.0],
+        [1.0, -0.7309547, -0.50616],
+        [-0.7309547, 1.0, -0.9287013],
+        [-0.50616, -0.9287013, 1.0],
     ]
 
     expected = pd.DataFrame(expected, columns=[0, 1, 2], index=[0, 1, 2])