Adding geomin to rotations

.coveragerc

@@ -1,4 +1,5 @@
 [run]
 source = factor_analyzer
 omit =
-    *factor_analyzer/analyze.py
+    *factor_analyzer/analyze.py
+    *tests/*

README.rst

@@ -82,6 +82,8 @@ and ``Rotator``:
     (e) quartimin (oblique rotation)
     (f) quartimax (orthogonal rotation)
     (g) equamax (orthogonal rotation)
+    (h) geomin_obl (oblique rotation)
+    (i) geomin_ort (orthogonal rotation)
 
 In adddition, the package includes a ``confirmatory_factor_analyzer``
 module with a stand-alone ``ConfirmatoryFactorAnalyzer`` class. The

doc/index.rst

@@ -73,6 +73,8 @@ and ``Rotator``:
     (e) quartimin (oblique rotation)
     (f) quartimax (orthogonal rotation)
     (g) equamax (orthogonal rotation)
+    (h) geomin_obl (oblique rotation)
+    (i) geomin_ort (orthogonal rotation)
 
 In adddition, the package includes a ``confirmatory_factor_analyzer``
 module with a stand-alone ``ConfirmatoryFactorAnalyzer`` class. The

factor_analyzer/rotator.py

@@ -12,9 +12,9 @@
 
 from sklearn.base import BaseEstimator
 
-ORTHOGONAL_ROTATIONS = ['varimax', 'oblimax', 'quartimax', 'equamax']
+ORTHOGONAL_ROTATIONS = ['varimax', 'oblimax', 'quartimax', 'equamax', 'geomin_ort']
 
-OBLIQUE_ROTATIONS = ['promax', 'oblimin', 'quartimin']
+OBLIQUE_ROTATIONS = ['promax', 'oblimin', 'quartimin', 'geomin_obl']
 
 POSSIBLE_ROTATIONS = ORTHOGONAL_ROTATIONS + OBLIQUE_ROTATIONS
 
@@ -36,6 +36,8 @@ class Rotator(BaseEstimator):
             (e) quartimin (oblique rotation)
             (f) quartimax (orthogonal rotation)
             (g) equamax (orthogonal rotation)
+            (h) geomin_obl (oblique rotation)
+            (i) geomin_ort (orthogonal rotation)
         Defaults to 'varimax'.
     normalize : bool or None, optional
         Whether to perform Kaiser normalization
@@ -56,6 +58,10 @@ class Rotator(BaseEstimator):
         The gamma level for the oblimin objective.
         Ignored if the method is not 'oblimin'.
         Defaults to 0.
+    delta : float, optional
+        The delta level for geomin objectives.
+         Ignored if the method is not 'geomin_*'.
+        Defaults to 0.01.
     max_iter : int, optional
         The maximum number of iterations.
         Used for varimax and oblique rotations.
@@ -112,6 +118,7 @@ def __init__(self,
                  power=4,
                  kappa=0,
                  gamma=0,
+                 delta=0.01,
                  max_iter=500,
                  tol=1e-5):
 
@@ -120,6 +127,7 @@ def __init__(self,
         self.power = power
         self.kappa = kappa
         self.gamma = gamma
+        self.delta = delta
         self.max_iter = max_iter
         self.tol = tol
 
@@ -256,6 +264,35 @@ def _equamax_obj(self, loadings):
         criterion = f1 + f2
         return {'grad': gradient, 'criterion': criterion}
 
+    def _geomin_obj(self, loadings):
+        """
+        Geomin function objective
+
+        Parameters
+        ----------
+        loadings : array-like
+            The loading matrix
+
+        Returns
+        -------
+        gradient_dict : dict
+            A dictionary with
+            - grad : np.array
+                The gradient.
+            - criterion : float
+                The value of the criterion for the objective.
+        """
+        p, k = loadings.shape
+
+        loadings2 = loadings**2 + self.delta
+
+        pro = np.exp(np.log(loadings2).sum(1) / k)
+        rep = np.repeat(pro, k, axis=0).reshape(p, k)
+
+        gradient = (2 / k) * (loadings / loadings2) * rep
+        criterion = np.sum(pro)
+        return {'grad': gradient, 'criterion': criterion}
+
     def _oblique(self, loadings, method):
         """
         A generic function for performing
@@ -286,6 +323,8 @@ def _oblique(self, loadings, method):
             objective = self._oblimin_obj
         elif method == 'quartimin':
             objective = self._quartimin_obj
+        elif method == 'geomin_obl':
+            objective = self._geomin_obj
 
         # initialize the rotation matrix
         n_rows, n_cols = loadings.shape
@@ -370,6 +409,8 @@ def _orthogonal(self, loadings, method):
             objective = self._quartimax_obj
         elif method == 'equamax':
             objective = self._equamax_obj
+        elif method == 'geomin_ort':
+            objective = self._geomin_obj
 
         arr = loadings.copy()
 

factor_analyzer/test_utils.py

@@ -228,9 +228,6 @@ def check_close(data1, data2, rel_tol=0.0, abs_tol=0.1,
         data1 = normalize(data1, absolute)
         data2 = normalize(data2, absolute)
 
-    # print(data1)
-    # print(data2)
-
     err_msg = 'r - py: {} != {}'
     assert data1.shape == data2.shape, err_msg.format(data1.shape, data2.shape)
 

tests/expected/test02/communalities_ml_geomin_obl_3_test02.csv

@@ -0,0 +1,11 @@
+"","x"
+"sex",0.568262782113056
+"zygosity",0.00335344116024749
+"moed",0.502211155728153
+"faed",0.732278623062979
+"faminc",0.331437013920703
+"english",0.640087242903545
+"math",0.609288716925608
+"socsci",0.748529483037574
+"natsci",0.641681114156784
+"vocab",0.73417183828982

tests/expected/test02/communalities_ml_geomin_ort_3_test02.csv

@@ -0,0 +1,11 @@
+"","x"
+"sex",0.568262782113056
+"zygosity",0.00335344116024749
+"moed",0.502211155728153
+"faed",0.732278623062979
+"faminc",0.331437013920703
+"english",0.640087242903545
+"math",0.609288716925608
+"socsci",0.748529483037574
+"natsci",0.641681114156784
+"vocab",0.73417183828982

tests/expected/test02/communalities_uls_geomin_obl_2_test02.csv

@@ -0,0 +1,11 @@
+"","x"
+"sex",0.0200698806283079
+"zygosity",0.00369309370387438
+"moed",0.502379851910451
+"faed",0.726577338166964
+"faminc",0.332505189916626
+"english",0.544999778173725
+"math",0.557021996057969
+"socsci",0.748062337667605
+"natsci",0.632540585595578
+"vocab",0.678479276307396

tests/expected/test02/communalities_uls_geomin_obl_3_test02.csv

@@ -0,0 +1,11 @@
+"","x"
+"sex",0.588608854233434
+"zygosity",0.00382390039248182
+"moed",0.504572482849861
+"faed",0.728337671409498
+"faminc",0.331860342110816
+"english",0.662097154411232
+"math",0.619120926251366
+"socsci",0.731949279298031
+"natsci",0.64930067481466
+"vocab",0.71150248325263

tests/expected/test02/communalities_uls_geomin_ort_2_test02.csv

@@ -0,0 +1,11 @@
+"","x"
+"sex",0.0200698806283079
+"zygosity",0.00369309370387438
+"moed",0.502379851910451
+"faed",0.726577338166964
+"faminc",0.332505189916626
+"english",0.544999778173725
+"math",0.557021996057969
+"socsci",0.748062337667605
+"natsci",0.632540585595578
+"vocab",0.678479276307396

tests/expected/test02/communalities_uls_geomin_ort_3_test02.csv

@@ -0,0 +1,11 @@
+"","x"
+"sex",0.588608854233434
+"zygosity",0.00382390039248182
+"moed",0.504572482849861
+"faed",0.728337671409498
+"faminc",0.331860342110816
+"english",0.662097154411232
+"math",0.619120926251366
+"socsci",0.731949279298031
+"natsci",0.64930067481466
+"vocab",0.71150248325263

tests/expected/test02/loading_ml_geomin_obl_3_test02.csv

@@ -0,0 +1,11 @@
+"","ML1","ML2","ML3"
+"sex",-0.000762424180361988,0.0144863568248402,0.753740579425805
+"zygosity",0.00836032084508804,0.0540337519005978,0.0120290443892905
+"moed",-0.0308305913041046,0.717662489907825,-0.0103105330331608
+"faed",0.00560377563048698,0.853885860518788,-0.00497216003434654
+"faminc",0.088593271872631,0.541812818111674,0.018664059359627
+"english",0.801373392520266,-0.0115490767443039,0.20620667237673
+"math",0.671162101786714,0.0538150112759279,-0.293868399277289
+"socsci",0.863835654302348,-0.0133638666020436,-0.0388515465301251
+"natsci",0.754721630451217,-0.0361046222604616,-0.220858975489651
+"vocab",0.849713016725188,0.0387564984246252,0.10398066569673

tests/expected/test02/loading_ml_geomin_ort_3_test02.csv

@@ -0,0 +1,11 @@
+"","ML1","ML2","ML3"
+"sex",-0.0494435245262262,0.0220695346868252,0.751886281483064
+"zygosity",0.0234129270234461,0.0519555635521807,0.0116562976387327
+"moed",0.181241090645935,0.685052508647751,-0.0079999269695229
+"faed",0.257357971369009,0.816104146267198,-0.00427750194975829
+"faminc",0.246685565300333,0.519964970597924,0.0149551078582187
+"english",0.78236066409413,0.00900079406183916,0.167086239703197
+"math",0.706706289662603,0.0631095064385724,-0.325381026154269
+"socsci",0.861415163462676,0.00597278536382461,-0.0803614744226611
+"natsci",0.758544289583419,-0.0201714914258842,-0.256679933020209
+"vocab",0.852624955056397,0.0570312190299659,0.0628474376734453

tests/expected/test02/loading_uls_geomin_obl_2_test02.csv

@@ -0,0 +1,11 @@
+"","MR1","MR2"
+"sex",-0.148097688436596,0.0640668247222744
+"zygosity",0.00926914066317219,0.0568602293413421
+"moed",-0.0293471965986073,0.717481062078967
+"faed",0.00641761627269584,0.850347742339245
+"faminc",0.0836351297182976,0.544788769783947
+"english",0.7323675050852,0.0179947840118985
+"math",0.737081067880239,0.0278811783103989
+"socsci",0.869669808042782,-0.015540390866951
+"natsci",0.810400558343853,-0.052922992443324
+"vocab",0.802492366408111,0.0608741581280668

tests/expected/test02/loading_uls_geomin_obl_3_test02.csv

@@ -0,0 +1,11 @@
+"","MR1","MR2","MR3"
+"sex",-0.0037951276006728,0.0170471822973247,0.766723858407271
+"zygosity",0.0132831696383983,0.0549321897093055,0.0152984160530353
+"moed",-0.0311667473841201,0.719374965767141,-0.0102226742311569
+"faed",0.00588990745158584,0.851500409638379,-0.00561207268328726
+"faminc",0.088647594887529,0.542272376359976,0.0196435943219139
+"english",0.814300188612467,-0.0167383523508025,0.22695076307264
+"math",0.685185995299752,0.0523767271171973,-0.279701210933423
+"socsci",0.853901182186794,-0.0112468992751859,-0.0361343611400812
+"natsci",0.764761045544525,-0.0341272144530272,-0.205266255147408
+"vocab",0.833734995516388,0.0464511125256671,0.102564983413607

tests/expected/test02/loading_uls_geomin_ort_2_test02.csv

@@ -0,0 +1,11 @@
+"","MR1","MR2"
+"sex",-0.130830186146755,0.054259950144959
+"zygosity",0.0244142922091219,0.0552715394436468
+"moed",0.16194647306451,0.690036213592894
+"faed",0.233083284970722,0.819903817844926
+"faminc",0.228750352674428,0.529305583551424
+"english",0.736227398973548,0.0543766362499309
+"math",0.743570309327008,0.064143654045391
+"socsci",0.864414712726003,0.0289977178030043
+"natsci",0.795256350990232,-0.0100292700934536
+"vocab",0.81769273813291,0.0992504308638485

tests/expected/test02/loading_uls_geomin_ort_3_test02.csv

@@ -0,0 +1,11 @@
+"","MR1","MR2","MR3"
+"sex",-0.0532444447811588,0.0246865000622711,0.76496441298773
+"zygosity",0.028284819500953,0.0529754424978699,0.0147220883331733
+"moed",0.180542636341071,0.686961216875825,-0.00711735663995621
+"faed",0.255960408725063,0.814135701739938,-0.00406140790000246
+"faminc",0.246131131884449,0.520580343886359,0.0163874268323194
+"english",0.792091541006628,0.00421988696762836,0.186199725078558
+"math",0.719413935634326,0.0618900911260461,-0.312630953912028
+"socsci",0.851929824850015,0.00742472376838442,-0.07814469484652
+"natsci",0.768221962653517,-0.0182577062908545,-0.242495462625807
+"vocab",0.838860805745795,0.0636801337752757,0.0613134554242502

tests/generate_r_output.r

@@ -0,0 +1,87 @@
+# install packages, if they're not already installed
+packages <- c('argparse', 'psych', 'GPArotation')
+new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
+if(length(new_packages)) {
+  install.packages(new_packages)
+}
+
+library('argparse')
+library('psych')
+
+# mapping list for fit methods
+mapping1 <- list(minres = 'uls',
+                 mle = 'ml')
+
+# mapping list for rotations
+mapping2 <- list(geominT = 'geomin_ort',
+                 geominQ = 'geomin_obl')
+
+# argument parser
+parser <- ArgumentParser(description='Fit some factor models')
+
+parser$add_argument('-n', '--n_factors', type='integer', nargs='+',
+                    default=2, help='integer(s) specifying the number of factors')
+
+parser$add_argument('-f', '--fit_methods', type='character', nargs='+',
+                    default='minres', help='Fit method(s)')
+
+parser$add_argument('-r', '--rotations', type="character", nargs='+',
+                    default='promax', help='Rotaton(s)')
+
+parser$add_argument('-t', '--test_file', type="character", nargs='+',
+                    default='test02.csv', help='Test file')
+
+
+
+# parse the arguments into a list
+args <- parser$parse_args()
+
+# get the input path and directory
+path <- args$test_file
+dir <- dirname(path)
+filename <- basename(path)
+
+# read in the data
+df <- read.csv(path)
+
+# loop through all of the conditions
+for (n in args$n_factors) {
+    for (fm in args$fit_methods) {
+        for (rot in args$rotations) {
+
+          # fit the factor model
+          res <- fa(cor(df),
+                    nfactors = n,
+                    n.obs = nrow(df),
+                    fm = fm,
+                    rotate = rot)
+
+          # get the correct method and rotation names
+          fm_name <- mapping1[[fm]]
+          fm_name <- if (length(fm_name) == 0) fm else fm_name
+
+          rot_name <- mapping2[[rot]]
+          rot_name <- if (length(rot_name) == 0) rot else rot_name
+
+          # write out the loadings
+          loadings_file <- paste('loading',
+                                 fm_name,
+                                 rot_name,
+                                 as.character(n),
+                                 filename,
+                                 sep='_')
+          loadings_file <- file.path(dir, loadings_file)
+          write.csv(res$loadings, loadings_file)
+
+          # write out the communalities
+          communalities_file <- paste('communalities',
+                                      fm_name,
+                                      rot_name,
+                                      as.character(n),
+                                      filename,
+                                      sep='_')
+          communalities_file <- file.path(dir, communalities_file)
+          write.csv(res$communalities, communalities_file)
+        }
+    }
+}