Repository: EducationalTestingService/factor_analyzer (language: python)

Issue:
Add pre-commit checks and apply them to all existing code
We should add the basic pre-commit checks (`black`, `pydocstyle`, `pylint`, etc.) since that would improve the quality and consistency of any new commits being added. We should also apply them to all of the existing code and reformat/fix any identified issues. 

We should also clean up the docstrings to update dates and email addresses and the like. 


Discussion / hints:

Fix patch:
diff --git a/.coveragerc b/.coveragerc
index 5b2a135..f373724 100644
--- a/.coveragerc
+++ b/.coveragerc
@@ -2,4 +2,4 @@
 source = factor_analyzer
 omit =
     *factor_analyzer/analyze.py
-    *tests/*
\ No newline at end of file
+    *tests/*
diff --git a/.gitignore b/.gitignore
index 8a6f03c..a3f68fe 100644
--- a/.gitignore
+++ b/.gitignore
@@ -11,4 +11,4 @@
 *static/
 *htmlcov/
 *templates/
-__pycache__
\ No newline at end of file
+__pycache__
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 0000000..46c507a
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,42 @@
+# See https://pre-commit.com for more information
+# See https://pre-commit.com/hooks.html for more git_hooks
+
+repos:
+-   repo: https://github.com/pre-commit/pre-commit-hooks
+    rev: v4.1.0
+    hooks:
+    -   id: trailing-whitespace
+    -   id: end-of-file-fixer
+        exclude: .(tsv|csv|json|json)$
+    -   id: check-added-large-files
+        args: [--maxkb=5000]
+    -   id: check-ast
+    -   id: check-json
+    -   id: debug-statements
+    -   id: check-merge-conflict
+-   repo: https://gitlab.com/pycqa/flake8
+    rev: 4.0.1
+    hooks:
+    - id: flake8
+      args: [--max-line-length=101]
+      exclude: ^(factor_analyzer/__init__.py)
+-   repo: https://github.com/PyCQA/isort
+    rev: 5.10.1
+    hooks:
+    -   id: isort
+        args: ["--profile", "black", "--filter-files"]
+-   repo: https://github.com/ikamensh/flynt/
+    rev: '0.76'
+    hooks:
+    -   id: flynt
+-   repo: https://github.com/psf/black
+    rev: 22.3.0
+    hooks:
+    -   id: black
+-   repo: https://github.com/pycqa/pydocstyle
+    rev: 6.1.1
+    hooks:
+    -   id: pydocstyle
+        args:
+        - --convention=pep257
+        additional_dependencies: ["toml"]
diff --git a/MANIFEST.in b/MANIFEST.in
index f439d04..a719767 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,4 +1,4 @@
 include LICENSE
 include README
 include *.rst
-include requirements.txt
\ No newline at end of file
+include requirements.txt
diff --git a/README.rst b/README.rst
index c1f942f..38a8515 100644
--- a/README.rst
+++ b/README.rst
@@ -2,27 +2,33 @@ FactorAnalyzer
 --------------
 
 .. image:: https://gitlab.com/EducationalTestingService/factor_analyzer/badges/main/pipeline.svg
-   :alt: Build status
    :target: https://gitlab.com/EducationalTestingService/factor_analyzer/-/pipelines
+   :alt: Build status
 
 .. image:: https://codecov.io/gh/EducationalTestingService/factor_analyzer/branch/main/graph/badge.svg
-  :target: https://codecov.io/gh/EducationalTestingService/factor_analyzer
+   :target: https://codecov.io/gh/EducationalTestingService/factor_analyzer
+   :alt: Code coverage
 
 .. image:: https://anaconda.org/ets/factor_analyzer/badges/installer/conda.svg
-    :target: https://anaconda.org/ets/factor_analyzer/
+   :target: https://anaconda.org/ets/factor_analyzer/
+   :alt: Conda version
 
 .. image:: https://img.shields.io/readthedocs/factor_analyzer/latest.svg
    :target: https://factor-analyzer.readthedocs.io/
    :alt: Docs
 
+.. image:: https://img.shields.io/badge/pre--commit-enabled-brightgreen?logo=pre-commit&logoColor=white
+   :target: https://github.com/pre-commit/pre-commit
+   :alt: Pre-commit checks
+
 
-This is a Python module to perform exploratory and factor analysis (EFA), with several 
+This is a Python module to perform exploratory and factor analysis (EFA), with several
 optional rotations. It also includes a class to perform confirmatory factor
-analysis (CFA), with certain pre-defined constraints. In expoloratory factor analysis,
+analysis (CFA), with certain pre-defined constraints. In exploratory factor analysis,
 factor extraction can be performed using a variety of estimation techniques. The
-``factor_analyzer`` package allows users to perfrom EFA using either (1) a minimum
+``factor_analyzer`` package allows users to perform EFA using either (1) a minimum
 residual (MINRES) solution, (2) a maximum likelihood (ML) solution, or (3) a principal
-factor solution. However, CFA can only be performe using an ML solution.
+factor solution. However, CFA can only be performed using an ML solution.
 
 Both the EFA and CFA classes within this package are fully compatible with `scikit-learn`.
 Portions of this code are ported from the excellent R library `psych`, and the `sem`
@@ -71,12 +77,12 @@ Two common types of rotations are:
    correlated.
 
 This package includes a ``factor_analyzer`` module with a stand-alone
-``FactorAnalyzer`` class. The class includes ``fit()`` and ``transform()`` 
+``FactorAnalyzer`` class. The class includes ``fit()`` and ``transform()``
 methods that enable users to perform factor analysis and score new data
-using the fitted factor model. Users can also perform optional otations
+using the fitted factor model. Users can also perform optional rotations
 on a factor loading matrix using the ``Rotator`` class.
 
-The following rotation options are available in both ``FactorAnalyzer`` 
+The following rotation options are available in both ``FactorAnalyzer``
 and ``Rotator``:
 
     (a) varimax (orthogonal rotation)
@@ -89,7 +95,7 @@ and ``Rotator``:
     (h) geomin_obl (oblique rotation)
     (i) geomin_ort (orthogonal rotation)
 
-In adddition, the package includes a ``confirmatory_factor_analyzer``
+In addition, the package includes a ``confirmatory_factor_analyzer``
 module with a stand-alone ``ConfirmatoryFactorAnalyzer`` class. The
 class includes ``fit()`` and ``transform()``  that enable users to perform
 confirmatory factor analysis and score new data using the fitted model.
@@ -107,21 +113,21 @@ Exploratory factor analysis example.
 
 .. code:: python
 
-  In [1]: import pandas as pd 
-     ...: from factor_analyzer import FactorAnalyzer                                                                                                     
+  In [1]: import pandas as pd
+     ...: from factor_analyzer import FactorAnalyzer
 
-  In [2]: df_features = pd.read_csv('tests/data/test02.csv')                                                                                             
+  In [2]: df_features = pd.read_csv('tests/data/test02.csv')
 
-  In [3]: fa = FactorAnalyzer(rotation=None)                                                                                                             
+  In [3]: fa = FactorAnalyzer(rotation=None)
 
-  In [4]: fa.fit(df_features)                                                                                                                            
-  Out[4]: 
+  In [4]: fa.fit(df_features)
+  Out[4]:
   FactorAnalyzer(bounds=(0.005, 1), impute='median', is_corr_matrix=False,
                  method='minres', n_factors=3, rotation=None, rotation_kwargs={},
                  use_smc=True)
 
-  In [5]: fa.loadings_                                                                                                                                   
-  Out[5]: 
+  In [5]: fa.loadings_
+  Out[5]:
   array([[-0.12991218,  0.16398151,  0.73823491],
          [ 0.03899558,  0.04658425,  0.01150343],
          [ 0.34874135,  0.61452341, -0.07255666],
@@ -133,8 +139,8 @@ Exploratory factor analysis example.
          [ 0.76041819, -0.23768727, -0.12068582],
          [ 0.81533404, -0.12494695,  0.17639684]])
 
-  In [6]: fa.get_communalities()                                                                                                                         
-  Out[6]: 
+  In [6]: fa.get_communalities()
+  Out[6]:
   array([0.5887579 , 0.00382308, 0.50452402, 0.72841182, 0.33184336,
          0.66208429, 0.61911037, 0.73194557, 0.64929612, 0.71149718])
 
@@ -142,24 +148,24 @@ Confirmatory factor analysis example.
 
 .. code:: python
 
-  In [1]: import pandas as pd                                                                                                                            
+  In [1]: import pandas as pd
 
-  In [2]: from factor_analyzer import (ConfirmatoryFactorAnalyzer, 
-     ...:                              ModelSpecificationParser)                                                                                         
+  In [2]: from factor_analyzer import (ConfirmatoryFactorAnalyzer,
+     ...:                              ModelSpecificationParser)
 
-  In [3]: df_features = pd.read_csv('tests/data/test11.csv')                                                                                             
+  In [3]: df_features = pd.read_csv('tests/data/test11.csv')
 
-  In [4]: model_dict = {"F1": ["V1", "V2", "V3", "V4"], 
-     ...:               "F2": ["V5", "V6", "V7", "V8"]} 
+  In [4]: model_dict = {"F1": ["V1", "V2", "V3", "V4"],
+     ...:               "F2": ["V5", "V6", "V7", "V8"]}
   In [5]: model_spec = ModelSpecificationParser.parse_model_specification_from_dict(df_features,
      ...:                                                                           model_dict)
 
-  In [6]: cfa = ConfirmatoryFactorAnalyzer(model_spec, disp=False)                                                                                       
+  In [6]: cfa = ConfirmatoryFactorAnalyzer(model_spec, disp=False)
 
-  In [7]: cfa.fit(df_features.values)                                                                                                                    
+  In [7]: cfa.fit(df_features.values)
 
-  In [8]: cfa.loadings_                                                                                                                                  
-  Out[8]: 
+  In [8]: cfa.loadings_
+  Out[8]:
   array([[0.99131285, 0.        ],
          [0.46074919, 0.        ],
          [0.3502267 , 0.        ],
@@ -169,13 +175,13 @@ Confirmatory factor analysis example.
          [0.        , 0.37602988],
          [0.        , 0.50049507]])
 
-  In [9]: cfa.factor_varcovs_                                                                                                                           
-  Out[9]: 
+  In [9]: cfa.factor_varcovs_
+  Out[9]:
   array([[1.        , 0.17385704],
          [0.17385704, 1.        ]])
 
-  In [10]: cfa.transform(df_features.values)                                                                                                             
-  Out[10]: 
+  In [10]: cfa.transform(df_features.values)
+  Out[10]:
   array([[-0.46852166, -1.08708035],
          [ 2.59025301,  1.20227783],
          [-0.47215977,  2.65697245],
@@ -192,12 +198,13 @@ Requirements
 -  ``pandas``
 -  ``scipy``
 -  ``scikit-learn``
+-  ``pre-commit``
 
 Contributing
 ------------
 
 Contributions to ``factor_analyzer`` are very welcome. Please file an issue
-on GitHub, or contact jbiggs@ets.org if you would like to contribute.
+in the repository if you would like to contribute.
 
 Installation
 ------------
diff --git a/doc/Makefile b/doc/Makefile
index aa26098..63a0a00 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -17,4 +17,4 @@ help:
 # Catch-all target: route all unknown targets to Sphinx using the new
 # "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
 %: Makefile
-	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
\ No newline at end of file
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
diff --git a/doc/conf.py b/doc/conf.py
index ec806a1..898c6ac 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
+"""Configure sphinx documentation."""
 # factor_analyzer documentation build configuration file, created by
 # sphinx-quickstart on Wed Dec  6 21:35:47 2017.
 #
@@ -19,7 +20,10 @@
 #
 import os
 import sys
-sys.path.insert(0, os.path.abspath('..'))
+
+import sphinx_rtd_theme
+
+sys.path.insert(0, os.path.abspath(".."))
 
 
 # -- General configuration ------------------------------------------------
@@ -31,59 +35,61 @@
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
-extensions = ['sphinx.ext.autodoc',
-              'sphinx.ext.viewcode',
-              'sphinx.ext.mathjax',
-              'sphinx.ext.napoleon',
-              'sphinx.ext.doctest',
-              'sphinx.ext.intersphinx',
-              'sphinx.ext.todo',
-              'sphinx.ext.coverage',
-              'sphinx.ext.mathjax',
-              'sphinx.ext.ifconfig',
-              'sphinx.ext.viewcode',
-              'sphinx.ext.githubpages']
+extensions = [
+    "sphinx.ext.autodoc",
+    "sphinx.ext.viewcode",
+    "sphinx.ext.mathjax",
+    "sphinx.ext.napoleon",
+    "sphinx.ext.doctest",
+    "sphinx.ext.intersphinx",
+    "sphinx.ext.todo",
+    "sphinx.ext.coverage",
+    "sphinx.ext.mathjax",
+    "sphinx.ext.ifconfig",
+    "sphinx.ext.viewcode",
+    "sphinx.ext.githubpages",
+]
 
 # Add any paths that contain templates here, relative to this directory.
-templates_path = ['_templates']
+templates_path = ["_templates"]
 
 # The suffix(es) of source filenames.
 # You can specify multiple suffix as a list of string:
 #
 # source_suffix = ['.rst', '.md']
-source_suffix = '.rst'
+source_suffix = ".rst"
 
 # The master toctree document.
-master_doc = 'index'
+master_doc = "index"
 
 # General information about the project.
-project = 'factor_analyzer'
-copyright = '2019, Jeremy Biggs'
-author = 'Jeremy Biggs'
+project = "factor_analyzer"
+copyright = "2019, Jeremy Biggs"
+author = "Jeremy Biggs"
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
 # built documents.
 #
 # The short X.Y version.
-version = '0.3.1'
+version = "0.4.0"
 # The full version, including alpha/beta/rc tags.
-release = '0.3.1'
+release = "0.4.0"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
 #
 # This is also used if you do content translation via gettext catalogs.
 # Usually you set "language" from the command line for these cases.
-language = None
+language = "en"
 
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
 # This patterns also effect to html_static_path and html_extra_path
-exclude_patterns = ['Thumbs.db', '.DS_Store']
+exclude_patterns = ["Thumbs.db", ".DS_Store"]
 
 # The name of the Pygments (syntax highlighting) style to use.
-pygments_style = 'sphinx'
+pygments_style = "sphinx"
 
 # If true, `todo` and `todoList` produce output, else they produce nothing.
 todo_include_todos = True
@@ -107,8 +113,8 @@
 # so a file named "default.css" will overwrite the builtin "default.css".
 # html_static_path = ['_static']
 
-import sphinx_rtd_theme
-html_theme = 'sphinx_rtd_theme'
+
+html_theme = "sphinx_rtd_theme"
 html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 
 # Custom sidebar templates, must be a dictionary that maps document names
@@ -117,12 +123,12 @@
 # This is required for the alabaster theme
 # refs: http://alabaster.readthedocs.io/en/latest/installation.html#sidebars
 html_sidebars = {
-    '**': [
-        'about.html',
-        'navigation.html',
-        'relations.html',  # needs 'show_related': True theme option to display
-        'searchbox.html',
-        'donate.html',
+    "**": [
+        "about.html",
+        "navigation.html",
+        "relations.html",  # needs 'show_related': True theme option to display
+        "searchbox.html",
+        "donate.html",
     ]
 }
 
@@ -130,7 +136,7 @@
 # -- Options for HTMLHelp output ------------------------------------------
 
 # Output file base name for HTML help builder.
-htmlhelp_basename = 'factor_analyzerdoc'
+htmlhelp_basename = "factor_analyzerdoc"
 
 
 # -- Options for LaTeX output ---------------------------------------------
@@ -139,15 +145,12 @@
     # The paper size ('letterpaper' or 'a4paper').
     #
     # 'papersize': 'letterpaper',
-
     # The font size ('10pt', '11pt' or '12pt').
     #
     # 'pointsize': '10pt',
-
     # Additional stuff for the LaTeX preamble.
     #
     # 'preamble': '',
-
     # Latex figure (float) alignment
     #
     # 'figure_align': 'htbp',
@@ -157,8 +160,13 @@
 # (source start file, target name, title,
 #  author, documentclass [howto, manual, or own class]).
 latex_documents = [
-    (master_doc, 'factor_analyzer.tex', 'factor\\_analyzer Documentation',
-     'Jeremy Biggs', 'manual'),
+    (
+        master_doc,
+        "factor_analyzer.tex",
+        "factor\\_analyzer Documentation",
+        "Jeremy Biggs",
+        "manual",
+    ),
 ]
 
 
@@ -167,8 +175,7 @@
 # One entry per manual page. List of tuples
 # (source start file, name, description, authors, manual section).
 man_pages = [
-    (master_doc, 'factor_analyzer', 'factor_analyzer Documentation',
-     [author], 1)
+    (master_doc, "factor_analyzer", "factor_analyzer Documentation", [author], 1)
 ]
 
 
@@ -178,9 +185,15 @@
 # (source start file, target name, title, author,
 #  dir menu entry, description, category)
 texinfo_documents = [
-    (master_doc, 'factor_analyzer', 'factor_analyzer Documentation',
-     author, 'factor_analyzer', 'One line description of project.',
-     'Miscellaneous'),
+    (
+        master_doc,
+        "factor_analyzer",
+        "factor_analyzer Documentation",
+        author,
+        "factor_analyzer",
+        "One line description of project.",
+        "Miscellaneous",
+    ),
 ]
 
 # -- Options for Epub output ----------------------------------------------
@@ -201,7 +214,7 @@
 # epub_uid = ''
 
 # A list of files that should not be packed into the epub file.
-epub_exclude_files = ['search.html']
+epub_exclude_files = ["search.html"]
 
 # Example configuration for intersphinx: refer to the Python standard library.
-intersphinx_mapping = {'https://docs.python.org/': None}
+intersphinx_mapping = {"https://docs.python.org/": None}
diff --git a/doc/factor_analyzer.rst b/doc/factor_analyzer.rst
index b3b92d1..92307ba 100644
--- a/doc/factor_analyzer.rst
+++ b/doc/factor_analyzer.rst
@@ -1,30 +1,32 @@
+.. _factor_analyzer_api:
+
 factor\_analyzer package
 ========================
 
-:py:mod:`factor\_analyzer\.analyze` Module
--------------------------------------------------
-.. automodule:: factor_analyzer.analyze
-    :members:
-    :undoc-members:
-    :show-inheritance:
-
 :py:mod:`factor\_analyzer\.factor\_analyzer` Module
---------------------------------------------------------
+---------------------------------------------------
 .. automodule:: factor_analyzer.factor_analyzer
     :members:
     :undoc-members:
     :show-inheritance:
 
 :py:mod:`factor\_analyzer\.confirmatory\_factor\_analyzer` Module
---------------------------------------------------------
+-----------------------------------------------------------------
 .. automodule:: factor_analyzer.confirmatory_factor_analyzer
     :members:
     :undoc-members:
     :show-inheritance:
 
 :py:mod:`factor\_analyzer\.rotator` Module
---------------------------------------------------------
+------------------------------------------
 .. automodule:: factor_analyzer.rotator
     :members:
     :undoc-members:
-    :show-inheritance:
\ No newline at end of file
+    :show-inheritance:
+
+:py:mod:`factor\_analyzer\.utils` Module
+----------------------------------------
+.. automodule:: factor_analyzer.utils
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/doc/index.rst b/doc/index.rst
index ef416f2..aa07880 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -6,17 +6,18 @@
 Welcome to the FactorAnalyzer documentation!
 ==============================================
 
-This is a Python module to perform exploratory and factor analysis (EFA), with several 
-optional rotations. It also includes a class to perform confirmatory factor
-analysis (CFA), with certain pre-defined constraints. In expoloratory factor analysis,
-factor extraction can be performed using a variety of estimation techniques. The
-``factor_analyzer`` package allows users to perfrom EFA using either (1) a minimum
-residual (MINRES) solution, (2) a maximum likelihood (ML) solution, or (3) a principal
-factor solution. However, CFA can only be performe using an ML solution.
-
-Both the EFA and CFA classes within this package are fully compatible with `scikit-learn`.
-Portions of this code are ported from the excellent R library `psych`, and the `sem`
-package provided inspiration for the CFA class.
+This is a Python module to perform exploratory and factor analysis (EFA), with
+several optional rotations. It also includes a class to perform confirmatory
+factor analysis (CFA), with certain pre-defined constraints. In exploratory
+factor analysis, factor extraction can be performed using a variety of
+estimation techniques. The :ref:`factor_analyzer <factor_analyzer_api>` package
+allows users to perform EFA using either (1) a minimum residual (MINRES)
+solution, (2) a maximum likelihood (ML) solution, or (3) a principal factor
+solution. However, CFA can only be performed using an ML solution.
+
+Both the EFA and CFA classes within this package are fully compatible with
+``scikit-learn``. Portions of this code are ported from the excellent R library
+``psych``, and the ``sem`` package provided inspiration for the CFA class.
 
 Description
 ============
@@ -30,7 +31,7 @@ from an EFA model describes the underlying relationships between each
 variable and the latent factors.
 
 Confirmatory factor analysis (CFA), a closely associated technique, is
-used to test an a priori hypothesis about latent relationships among sets
+used to test a priori hypothesis about latent relationships among sets
 of observed variables. In CFA, the researcher specifies the expected pattern
 of factor loadings (and possibly other constraints), and fits a model according
 to this specification.
@@ -57,14 +58,17 @@ Two common types of rotations are:
    upon the varimax rotation, but ultimately allows factors to become
    correlated.
 
-This package includes a ``factor_analyzer`` module with a stand-alone
-``FactorAnalyzer`` class. The class includes ``fit()`` and ``transform()`` 
-methods that enable users to perform factor analysis and score new data
-using the fitted factor model. Users can also perform optional rotations
-on a factor loading matrix using the ``Rotator`` class.
+This package includes a :mod:`~factor_analyzer.factor_analyzer` module with a
+stand-alone :class:`~factor_analyzer.factor_analyzer.FactorAnalyzer` class. The
+class includes :meth:`~factor_analyzer.factor_analyzer.FactorAnalyzer.fit()`
+and :meth:`~factor_analyzer.factor_analyzer.FactorAnalyzer.transform()`
+methods that enable users to perform factor analysis and score new data using
+the fitted factor model. Users can also perform optional rotations on a factor
+loading matrix using the :class:`~factor_analyzer.rotator.Rotator` class.
 
-The following rotation options are available in both ``FactorAnalyzer`` 
-and ``Rotator``:
+The following rotation options are available in both
+:class:`~factor_analyzer.factor_analyzer.FactorAnalyzer` and
+:class:`~factor_analyzer.rotator.Rotator`:
 
     (a) varimax (orthogonal rotation)
     (b) promax (oblique rotation)
@@ -76,21 +80,28 @@ and ``Rotator``:
     (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
-class includes ``fit()`` and ``transform()``  that enable users to perform
-confirmatory factor analysis and score new data using the fitted model.
-Performing CFA requires users to specify in advance a model specification
-with the expected factor loading relationships. This can be done using
-the ``ModelSpecificationParser`` class.
+In addition, the package includes a
+:mod:`~factor_analyzer.confirmatory_factor_analyzer` module with a stand-alone
+:class:`~factor_analyzer.confirmatory_factor_analyzer.ConfirmatoryFactorAnalyzer`
+class. The class includes
+:meth:`~factor_analyzer.confirmatory_factor_analyzer.ConfirmatoryFactorAnalyzer.fit()`
+and
+:meth:`~factor_analyzer.confirmatory_factor_analyzer.ConfirmatoryFactorAnalyzer.transform()`
+that enable users to perform confirmatory factor analysis and score new data
+using the fitted model. Performing CFA requires users to specify in advance a
+model specification with the expected factor loading relationships. This can be
+done using the
+:class:`~factor_analyzer.confirmatory_factor_analyzer.ModelSpecificationParser`
+class.
 
 Requirements
 ==================
--  Python 3.4 or higher
+-  Python 3.7 or higher
 -  ``numpy``
 -  ``pandas``
 -  ``scipy``
 -  ``scikit-learn``
+-  ``pre-commit``
 
 Installation
 ==================
@@ -109,10 +120,6 @@ Alternatively, you can install via ``conda`` with:
 
    factor_analyzer
 
-   confirmatory_factor_analyzer
-
-   rotator
-
 
 Indices and tables
 ==================
diff --git a/factor_analyzer/__init__.py b/factor_analyzer/__init__.py
index 03028ba..2ecab37 100644
--- a/factor_analyzer/__init__.py
+++ b/factor_analyzer/__init__.py
@@ -1,24 +1,35 @@
-# License: GPL2
 """
-:author: Jeremy Biggs (jbiggs@ets.org)
+This module performs exploratory and confirmatory factor analyses.
+
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
 :author: Nitin Madnani (nmadnani@ets.org)
 :organization: Educational Testing Service
-:date: 2021-10-18
+:date: 2022-09-05
 """
 
-from .confirmatory_factor_analyzer import ConfirmatoryFactorAnalyzer, ModelSpecification, ModelSpecificationParser
-from .factor_analyzer import FactorAnalyzer, calculate_bartlett_sphericity, calculate_kmo
+from .confirmatory_factor_analyzer import (
+    ConfirmatoryFactorAnalyzer,
+    ModelSpecification,
+    ModelSpecificationParser,
+)
+from .factor_analyzer import (
+    FactorAnalyzer,
+    calculate_bartlett_sphericity,
+    calculate_kmo,
+)
 from .rotator import Rotator
-from .utils import (commutation_matrix,
-                    corr,
-                    cov,
-                    covariance_to_correlation,
-                    duplication_matrix,
-                    duplication_matrix_pre_post,
-                    fill_lower_diag,
-                    get_symmetric_lower_idxs,
-                    get_symmetric_upper_idxs,
-                    impute_values,
-                    merge_variance_covariance,
-                    partial_correlations,
-                    smc)
+from .utils import (
+    commutation_matrix,
+    corr,
+    cov,
+    covariance_to_correlation,
+    duplication_matrix,
+    duplication_matrix_pre_post,
+    fill_lower_diag,
+    get_symmetric_lower_idxs,
+    get_symmetric_upper_idxs,
+    impute_values,
+    merge_variance_covariance,
+    partial_correlations,
+    smc,
+)
diff --git a/factor_analyzer/confirmatory_factor_analyzer.py b/factor_analyzer/confirmatory_factor_analyzer.py
index 85dcd9a..8dd93c0 100644
--- a/factor_analyzer/confirmatory_factor_analyzer.py
+++ b/factor_analyzer/confirmatory_factor_analyzer.py
@@ -1,10 +1,10 @@
 """
 Confirmatory factor analysis using machine learning methods.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
 :author: Nitin Madnani (nmadnani@ets.org)
 :organization: Educational Testing Service
-:date: 2021-10-18
+:date: 2022-09-05
 """
 
 import warnings
@@ -18,74 +18,47 @@
 from sklearn.utils import check_array
 from sklearn.utils.validation import check_is_fitted
 
-from .utils import (commutation_matrix,
-                    cov,
-                    covariance_to_correlation,
-                    duplication_matrix_pre_post,
-                    get_free_parameter_idxs,
-                    get_symmetric_lower_idxs,
-                    get_symmetric_upper_idxs,
-                    impute_values,
-                    merge_variance_covariance,
-                    unique_elements)
+from .utils import (
+    commutation_matrix,
+    cov,
+    covariance_to_correlation,
+    duplication_matrix_pre_post,
+    get_free_parameter_idxs,
+    get_symmetric_lower_idxs,
+    get_symmetric_upper_idxs,
+    impute_values,
+    merge_variance_covariance,
+    unique_elements,
+)
 
 
 class ModelSpecification:
     """
-    A class to encapsulate the model specification
-    for CFA. This class contains a number of specification
-    properties that are used in the CFA procedure.
+    Encapsulate the model specification for CFA.
+
+    This class contains a number of specification properties
+    that are used in the CFA procedure.
 
     Parameters
     ----------
     loadings : array-like
         The factor loadings specification.
-    error_vars : array-like
-        The error variance specification
-    factor_covs : array-like
-        The factor covariance specification.
-    factor_names : list of str or None
-        A list of factor names, if available.
-        Defaults to None.
-    variable_names : list of str or None
-        A list of variable names, if available.
-        Defaults to None.
-
-    Attributes
-    ----------
-    loadings : numpy array
-        The factor loadings specification.
-    error_vars : numpy array
-        The error variance specification
-    factor_covs : numpy array
-        The factor covariance specification.
     n_factors : int
         The number of factors.
     n_variables : int
         The number of variables.
-    n_lower_diag :int
-        The number of elements in the `factor_covs`
-        array, which is equal to the lower diagonal
-        of the factor covariance matrix.
-    loadings_free : numpy array
-        The indexes of "free" factor loading parameters.
-    error_vars_free : numpy array
-        The indexes of "free" error variance parameters.
-    factor_covs_free : numpy array
-        The indexes of "free" factor covariance parameters.
     factor_names : list of str or None
         A list of factor names, if available.
+        Defaults to ``None``.
     variable_names : list of str or None
         A list of variable names, if available.
+        Defaults to ``None``.
     """
 
-    def __init__(self,
-                 loadings,
-                 n_factors,
-                 n_variables,
-                 factor_names=None,
-                 variable_names=None):
-
+    def __init__(
+        self, loadings, n_factors, n_variables, factor_names=None, variable_names=None
+    ):
+        """Initialize the specification."""
         assert isinstance(loadings, np.ndarray)
         assert loadings.shape[0] == n_variables
         assert loadings.shape[1] == n_factors
@@ -107,53 +80,66 @@ def __init__(self,
         self._factor_covs_free = get_symmetric_lower_idxs(n_factors, False)
 
     def __str__(self):
-        return '<ModelSpecification object at {}>'.format(hex(id(self)))
+        """Represent model specification as a string."""
+        return f"<ModelSpecification object at {hex(id(self))}>"
 
     def copy(self):
+        """Return a copy of the model specification."""
         return deepcopy(self)
 
     @property
     def loadings(self):
+        """Get the factor loadings specification."""
         return self._loadings.copy()
 
     @property
     def error_vars(self):
+        """Get the error variance specification."""
         return self._error_vars.copy()
 
     @property
     def factor_covs(self):
+        """Get the factor covariance specification."""
         return self._factor_covs.copy()
 
     @property
     def loadings_free(self):
+        """Get the indices of "free" factor loading parameters."""
         return self._loadings_free.copy()
 
     @property
     def error_vars_free(self):
+        """Get the indices of "free" error variance parameters."""
         return self._error_vars_free.copy()
 
     @property
     def factor_covs_free(self):
+        """Get the indices of "free" factor covariance parameters."""
         return self._factor_covs_free.copy()
 
     @property
     def n_variables(self):
+        """Get the number of variables."""
         return self._n_variables
 
     @property
     def n_factors(self):
+        """Get the number of factors."""
         return self._n_factors
 
     @property
     def n_lower_diag(self):
+        """Get the lower diagonal of the factor covariance matrix."""
         return self._n_lower_diag
 
     @property
     def factor_names(self):
+        """Get list of factor names, if available."""
         return self._factor_names
 
     @property
     def variable_names(self):
+        """Get list of variable names, if available."""
         return self._variable_names
 
     def get_model_specification_as_dict(self):
@@ -166,56 +152,58 @@ def get_model_specification_as_dict(self):
             The model specification keys and values,
             as a dictionary.
         """
-        return {'loadings': self._loadings.copy(),
-                'error_vars': self._error_vars.copy(),
-                'factor_covs': self._factor_covs.copy(),
-                'loadings_free': self._loadings_free.copy(),
-                'error_vars_free': self._error_vars_free.copy(),
-                'factor_covs_free': self._factor_covs_free.copy(),
-                'n_variables': self._n_variables,
-                'n_factors': self._n_factors,
-                'n_lower_diag': self._n_lower_diag,
-                'variable_names': self._variable_names,
-                'factor_names': self._factor_names}
+        return {
+            "loadings": self._loadings.copy(),
+            "error_vars": self._error_vars.copy(),
+            "factor_covs": self._factor_covs.copy(),
+            "loadings_free": self._loadings_free.copy(),
+            "error_vars_free": self._error_vars_free.copy(),
+            "factor_covs_free": self._factor_covs_free.copy(),
+            "n_variables": self._n_variables,
+            "n_factors": self._n_factors,
+            "n_lower_diag": self._n_lower_diag,
+            "variable_names": self._variable_names,
+            "factor_names": self._factor_names,
+        }
 
 
 class ModelSpecificationParser:
     """
-    A class to generate the model specification for CFA.
-    This class includes two static methods to generate the
-    ``ModelSpecification`` object from either a dictionary
+    Generate the model specification for CFA.
+
+    This class includes two static methods to generate a
+    :class:`ModelSpecification` object from either a dictionary
     or a numpy array.
     """
 
     @staticmethod
     def parse_model_specification_from_dict(X, specification=None):
         """
-        Generate the model specification from a
-        dictionary. The keys in the dictionary
-        should be the factor names, and the values
-        should be the feature names. If this method
-        is used to create the ``ModelSpecification``,
-        then factor names and variable names will
-        be added as properties to that object.
+        Generate the model specification from a dictionary.
+
+        The keys in the dictionary should be the factor names, and the
+        values should be the feature names. If this method is used to
+        create the :class:`ModelSpecification` object, then factor names
+        and variable names will be added as properties to that object.
 
         Parameters
         ----------
         X : array-like
             The data set that will be used for CFA.
         specification : dict or None
-            A dictionary with the loading details. If None, the matrix will
+            A dictionary with the loading details. If ``None``, the matrix will
             be created assuming all variables load on all factors.
-            Defaults to None.
+            Defaults to ``None``.
 
         Returns
         -------
         ModelSpecification
-            A model specification object
+            A model specification object.
 
         Raises
         ------
         ValueError
-            If `specification` is not in the expected format.
+            If ``specification`` is not in the expected format.
 
         Examples
         --------
@@ -233,53 +221,63 @@ def parse_model_specification_from_dict(X, specification=None):
             loadings = np.ones((n_factors, n_factors), dtype=int)
         elif isinstance(specification, dict):
             factor_names = list(specification)
-            variable_names = unique_elements([v for f in specification.values() for v in f])
+            variable_names = unique_elements(
+                [v for f in specification.values() for v in f]
+            )
             loadings_new = {}
             for factor in factor_names:
-                loadings_for_factor = pd.Series(variable_names).isin(specification[factor])
+                loadings_for_factor = pd.Series(variable_names).isin(
+                    specification[factor]
+                )
                 loadings_for_factor = loadings_for_factor.astype(int)
                 loadings_new[factor] = loadings_for_factor
             loadings = pd.DataFrame(loadings_new).values
             n_variables, n_factors = loadings.shape
         else:
-            raise ValueError('The model `specification` must be either a dict '
-                             'or None, not {}'.format(type(specification)))
-
-        return ModelSpecification(**{'loadings': loadings,
-                                     'n_variables': n_variables,
-                                     'n_factors': n_factors,
-                                     'factor_names': factor_names,
-                                     'variable_names': variable_names})
+            raise ValueError(
+                "The model `specification` must be either a dict "
+                "or None, not {}".format(type(specification))
+            )
+
+        return ModelSpecification(
+            **{
+                "loadings": loadings,
+                "n_variables": n_variables,
+                "n_factors": n_factors,
+                "factor_names": factor_names,
+                "variable_names": variable_names,
+            }
+        )
 
     @staticmethod
     def parse_model_specification_from_array(X, specification=None):
         """
-        Generate the model specification from
-        an array. The columns should correspond to
-        the factors, and the rows should correspond to
-        the variables. If this method is used to create
-        the ``ModelSpecification``, then no factor names
-        and variable names will be added as properties
-        to that object.
+        Generate the model specification from a numpy array.
+
+        The columns should correspond to the factors, and the rows
+        should correspond to the variables. If this method is used
+        to create the :class:`ModelSpecification` object, then *no* factor
+        names and variable names will be added as properties to that
+        object.
 
         Parameters
         ----------
         X : array-like
             The data set that will be used for CFA.
         specification : array-like or None
-            An array with the loading details. If None, the matrix will
+            An array with the loading details. If ``None``, the matrix will
             be created assuming all variables load on all factors.
-            Defaults to None.
+            Defaults to ``None``.
 
         Returns
         -------
         ModelSpecification
-            A model specification object
+            A model specification object.
 
         Raises
         ------
         ValueError
-            If `specification` is not in the expected format.
+            If ``specification`` is not in the expected format.
 
         Examples
         --------
@@ -296,73 +294,70 @@ def parse_model_specification_from_array(X, specification=None):
             n_variables, n_factors = X.shape[1], X.shape[1]
             loadings = np.ones((n_factors, n_factors), dtype=int)
         elif isinstance(specification, (np.ndarray, pd.DataFrame)):
-                n_variables, n_factors = specification.shape
-                if isinstance(specification, pd.DataFrame):
-                    loadings = specification.values.copy()
-                else:
-                    loadings = specification.copy()
+            n_variables, n_factors = specification.shape
+            if isinstance(specification, pd.DataFrame):
+                loadings = specification.values.copy()
+            else:
+                loadings = specification.copy()
         else:
-            raise ValueError('The model `specification` must be either a numpy array '
-                             'or None, not {}'.format(type(specification)))
+            raise ValueError(
+                "The model `specification` must be either a numpy array "
+                "or None, not {}".format(type(specification))
+            )
 
-        return ModelSpecification(**{'loadings': loadings,
-                                     'n_variables': n_variables,
-                                     'n_factors': n_factors})
+        return ModelSpecification(
+            **{"loadings": loadings, "n_variables": n_variables, "n_factors": n_factors}
+        )
 
 
 class ConfirmatoryFactorAnalyzer(BaseEstimator, TransformerMixin):
     """
-    A ConfirmatoryFactorAnalyzer class, which fits a
-    confirmatory factor analysis model using maximum likelihood.
+    Fit a confirmatory factor analysis model using maximum likelihood.
 
     Parameters
     ----------
-    specification : ModelSpecificaition object or None, optional
-        A model specification. This must be a ``ModelSpecificaiton`` object
-        or None. If None, the ModelSpecification will be generated assuming
-        that n_factors == n_variables, and that all variables load on all
-        factors. Note that this could mean the factor model is not
-        identified, and the optimization could fail.
-        Defaults to None.
+    specification : ModelSpecification object or None, optional
+        A model specification. This must be a :class:`ModelSpecification` object
+        or ``None``. If ``None``, a :class:`ModelSpecification` object will be
+        generated assuming that ``n_factors`` == ``n_variables``, and that
+        all variables load on all factors. Note that this could mean the
+        factor model is not identified, and the optimization could fail.
+        Defaults to `None`.
     n_obs : int or None, optional
-        The number of observations in the original
-        data set. If this is not passed and `is_cov_matrix=True`,
+        The number of observations in the original data set.
+        If this is not passed and ``is_cov_matrix`` is ``True``,
         then an error will be raised.
-        Defaults to None.
+        Defaults to ``None``.
     is_cov_matrix : bool, optional
-        Whether the input `X` is a
-        covariance matrix. If False,
+        Whether the input ``X`` is a covariance matrix. If ``False``,
         assume it is the full data set.
-        Defaults to False.
+        Defaults to ``False``.
     bounds : list of tuples or None, optional
-        A list of minimum and maximum
-        boundaries for each element of the
-        input array. This must equal `x0`,
-        which is the input array from your
-        parsed and combined model specification.
+        A list of minimum and maximum boundaries for each element
+        of the input array. This must equal ``x0``, which is the
+        input array from your parsed and combined model specification.
+
         The length is:
-          ((n_factors * n_variables) + n_variables + n_factors + 
-          (((n_factors * n_factors) - n_factors) // 2)
-        If None, nothing will be bounded.
-        Defaults to None.
+        ((n_factors * n_variables) + n_variables + n_factors +
+        (((n_factors * n_factors) - n_factors) // 2)
+
+        If `None`, nothing will be bounded.
+        Defaults to ``None``.
     max_iter : int, optional
-        The maximum number of iterations
-        for the optimization routine.
+        The maximum number of iterations for the optimization routine.
         Defaults to 200.
     tol : float or None, optional
         The tolerance for convergence.
-        Defaults to None.
+        Defaults to ``None``.
     disp : bool, optional
-        Whether to print the scipy
-        optimization fmin message to
+        Whether to print the scipy optimization ``fmin`` message to
         standard output.
-        Defaults to True.
+        Defaults to ``True``.
 
     Raises
     ------
     ValueError
-        If ``is_cov_matrix`` is True, and n_obs
-        is not provided.
+        If `is_cov_matrix` is `True`, and `n_obs` is not provided.
 
     Attributes
     ----------
@@ -425,21 +420,25 @@ class ConfirmatoryFactorAnalyzer(BaseEstimator, TransformerMixin):
            [-0.27863554, -0.7695101 ]])
     """
 
-    def __init__(self,
-                 specification=None,
-                 n_obs=None,
-                 is_cov_matrix=False,
-                 bounds=None,
-                 max_iter=200,
-                 tol=None,
-                 impute='median',
-                 disp=True):
-
+    def __init__(
+        self,
+        specification=None,
+        n_obs=None,
+        is_cov_matrix=False,
+        bounds=None,
+        max_iter=200,
+        tol=None,
+        impute="median",
+        disp=True,
+    ):
+        """Initialize the analyzer object."""
         # if the input is going to be a covariance matrix, rather than
         # the full data set, then users must pass the number of observations
         if is_cov_matrix and n_obs is None:
-            raise ValueError('If `is_cov_matrix=True`, you must provide '
-                             'the number of observations, `n_obs`.')
+            raise ValueError(
+                "If `is_cov_matrix=True`, you must provide "
+                "the number of observations, `n_obs`."
+            )
 
         self.specification = specification
         self.n_obs = n_obs
@@ -465,57 +464,61 @@ def __init__(self,
         self._n_lower_diag = None
 
     @staticmethod
-    def _combine(loadings,
-                 error_vars,
-                 factor_vars,
-                 factor_covs,
-                 n_factors,
-                 n_variables,
-                 n_lower_diag):
+    def _combine(
+        loadings,
+        error_vars,
+        factor_vars,
+        factor_covs,
+        n_factors,
+        n_variables,
+        n_lower_diag,
+    ):
         """
-        Combine a set of multi-dimensional loading,
-        error variance, factor variance, and factor
-        covariance matrices into a one-dimensional
-        (X, 1) matrix, where the length of X is the
-        length of all elements across all of the matrices.
+        Combine given multi-dimensional matrics into a one-dimensional matrix.
+
+        Combine a set of multi-dimensional loading, error variance, factor
+        variance, and factor covariance matrices into a one-dimensional
+        (X, 1) matrix, where the length of X is the length of all elements
+        across all of the matrices.
 
         Parameters
         ----------
         loadings : array-like
-            The loadings matrix (n_factors * n_variables)
-        error_vars : array-like (n_variables * 1)
+            The loadings matrix (``n_factors`` * ``n_variables``)
+        error_vars : array-like (``n_variables`` * 1)
             The error variance array.
-        factor_vars : array-like (n_factors * 1)
+        factor_vars : array-like (``n_factors`` * 1)
             The factor variance array.
-        factor_covs : array-like (n_lower_diag * 1)
-            The factor covariance array
+        factor_covs : array-like (``n_lower_diag`` * 1)
+            The factor covariance array.
         n_factors : int
             The number of factors.
         n_variables : int
             The number of variables.
         n_lower_diag :int
-            The number of elements in the `factor_covs`
-            array, which is equal to the lower diagonal
-            of the factor covariance matrix.
+            The number of elements in the ``factor_covs`` array, which
+            is equal to the lower diagonal of the factor covariance
+            matrix.
 
         Returns
         -------
-        array : np.arrays
+        array : np.ndarray
             The combined (X, 1) array.
         """
-        loadings = loadings.reshape(n_factors * n_variables, 1, order='F')
-        error_vars = error_vars.reshape(n_variables, 1, order='F')
-        factor_vars = factor_vars.reshape(n_factors, 1, order='F')
-        factor_covs = factor_covs.reshape(n_lower_diag, 1, order='F')
+        loadings = loadings.reshape(n_factors * n_variables, 1, order="F")
+        error_vars = error_vars.reshape(n_variables, 1, order="F")
+        factor_vars = factor_vars.reshape(n_factors, 1, order="F")
+        factor_covs = factor_covs.reshape(n_lower_diag, 1, order="F")
         return np.concatenate([loadings, error_vars, factor_vars, factor_covs])
 
     @staticmethod
     def _split(x, n_factors, n_variables, n_lower_diag):
         """
-        Given a one-dimensional array, split it
-        into a set of multi-dimensional loading,
-        error variance, factor variance, and factor
-        covariance matrices.
+        Split given one-dimensional array into multi-dimensional arrays.
+
+        Given a one-dimensional array, split it into a set of
+        multi-dimensional loading, error variance, factor variance,
+        and factor covariance matrices.
 
         Parameters
         ----------
@@ -526,80 +529,91 @@ def _split(x, n_factors, n_variables, n_lower_diag):
         n_variables : int
             The number of variables.
         n_lower_diag : int
-            The number of elements in the `factor_covs`
-            array, which is equal to the lower.
+            The number of elements in the ``factor_covs`` array, which is
+            equal to the lower diagonal of the factor covariance matrix.
 
         Returns
         -------
         loadings : array-like
-            The loadings matrix (n_factors * n_variables)
-        error_vars : array-like (n_variables * 1)
+            The loadings matrix (``n_factors`` * ``n_variables``)
+        error_vars : array-like (``n_variables`` * 1)
             The error variance array.
-        factor_vars : array-like (n_factors * 1)
+        factor_vars : array-like (``n_factors`` * 1)
             The factor variance array.
-        factor_covs : array-like (n_lower_diag * 1)
-            The factor covariance array
+        factor_covs : array-like (``n_lower_diag`` * 1)
+            The factor covariance array.
         """
         loadings_ix = int(n_factors * n_variables)
         error_vars_ix = n_variables + loadings_ix
         factor_vars_ix = n_factors + error_vars_ix
         factor_covs_ix = n_lower_diag + factor_vars_ix
-        return (x[:loadings_ix].reshape((n_variables, n_factors), order='F'),
-                x[loadings_ix:error_vars_ix].reshape((n_variables, 1), order='F'),
-                x[error_vars_ix:factor_vars_ix].reshape((n_factors, 1), order='F'),
-                x[factor_vars_ix:factor_covs_ix].reshape((n_lower_diag, 1), order='F'))
+        return (
+            x[:loadings_ix].reshape((n_variables, n_factors), order="F"),
+            x[loadings_ix:error_vars_ix].reshape((n_variables, 1), order="F"),
+            x[error_vars_ix:factor_vars_ix].reshape((n_factors, 1), order="F"),
+            x[factor_vars_ix:factor_covs_ix].reshape((n_lower_diag, 1), order="F"),
+        )
 
     def _objective(self, x0, cov_mtx, loadings):
         """
-        The objective function.
+        Get the objective function for the analyzer.
 
         Parameters
         ----------
         x0: array-like
-            The combined (X, 1) array. These are the
-            initial values for the `minimize()` function.
+            The combined (X, 1) array. These are the initial values for
+            the ``minimize`()` function.
         cov_mtx : array-like
-            The covariance matrix from the original data
-            set.
+            The covariance matrix from the original data set.
         loadings : array-like
-            The loadings matrix (n_factors * n_variables)
-            from the model parser. This tells the objective
-            function what elements should be fixed.
+            The loadings matrix (``n_factors`` * ``n_variables``)
+            from the model parser. This tells the objective function what
+            elements should be fixed.
 
         Returns
         -------
         error : float
             The error from the objective function.
         """
-        (loadings_init,
-         error_vars_init,
-         factor_vars_init,
-         factor_covs_init) = self._split(x0,
-                                         self.model.n_factors,
-                                         self.model.n_variables,
-                                         self.model.n_lower_diag)
+        (
+            loadings_init,
+            error_vars_init,
+            factor_vars_init,
+            factor_covs_init,
+        ) = self._split(
+            x0, self.model.n_factors, self.model.n_variables, self.model.n_lower_diag
+        )
 
         # set the loadings to zero where applicable
         loadings_init[np.where(loadings == 0)] = 0
 
         # combine factor variances and covariances into a single matrix
-        factor_varcov_init = merge_variance_covariance(factor_vars_init, factor_covs_init)
+        factor_varcov_init = merge_variance_covariance(
+            factor_vars_init, factor_covs_init
+        )
 
         # make the error variance into a variance-covariance matrix
         error_varcov_init = merge_variance_covariance(error_vars_init)
 
         # make the factor variance-covariance matrix into a correlation matrix
-        with np.errstate(all='ignore'):
+        with np.errstate(all="ignore"):
             factor_varcov_init = covariance_to_correlation(factor_varcov_init)
 
         # calculate sigma-theta, needed for the objective function
-        sigma_theta = loadings_init.dot(factor_varcov_init) \
-                                   .dot(loadings_init.T) + error_varcov_init
-
-        with np.errstate(all='ignore'):
-            error = -(((-self.n_obs * self.model.n_variables / 2) * np.log(2 * np.pi)) -
-                      (self.n_obs / 2) * (np.log(np.linalg.det(sigma_theta)) +
-                                          np.trace(cov_mtx.dot(np.linalg.inv(sigma_theta)))))
+        sigma_theta = (
+            loadings_init.dot(factor_varcov_init).dot(loadings_init.T)
+            + error_varcov_init
+        )
+
+        with np.errstate(all="ignore"):
+            error = -(
+                ((-self.n_obs * self.model.n_variables / 2) * np.log(2 * np.pi))
+                - (self.n_obs / 2)
+                * (
+                    np.log(np.linalg.det(sigma_theta))
+                    + np.trace(cov_mtx.dot(np.linalg.inv(sigma_theta)))
+                )
+            )
 
             # make sure the error is greater than or
             # equal to zero before we return it; we
@@ -615,22 +629,18 @@ def fit(self, X, y=None):
         Parameters
         ----------
         X : array-like
-            The data to use for confirmatory
-            factor analysis. If this is just a
-            covariance matrix, make sure `is_cov_matrix`
-            was set to True.
+            The data to use for confirmatory factor analysis. If this is just a
+            covariance matrix, make sure ``is_cov_matrix`` was set to ``True``.
         y : ignored
 
         Raises
         ------
         ValueError
-            If the specification is not None or a
-            ``ModelSpecification`` object
+            If the specification is not None or a :class:`ModelSpecification` object.
         AssertionError
-            If ``is_cov_matrix=True`` and the matrix
-            is not square.
+            If ``is_cov_matrix`` was ``True`` and the matrix is not square.
         AssertionError
-            If len(bounds) != len(x0)
+            If ``len(bounds)`` != ``len(x0)``
 
         Examples
         --------
@@ -654,22 +664,23 @@ def fit(self, X, y=None):
                [0.        , 0.50049507]])
         """
         if self.specification is None:
-            self.model = ModelSpecificationParser.parse_model_specification_from_array(X)
+            self.model = ModelSpecificationParser.parse_model_specification_from_array(
+                X
+            )
         elif isinstance(self.specification, ModelSpecification):
             self.model = self.specification.copy()
         else:
-            raise ValueError('The `specification` must be None or `ModelSpecification` '
-                             'instance, not {}'.format(type(self.specification)))
+            raise ValueError(
+                "The `specification` must be None or `ModelSpecification` "
+                "instance, not {}".format(type(self.specification))
+            )
 
         if isinstance(X, pd.DataFrame):
             X = X.values
 
         # now check the array, and make sure it
         # meets all of our expected criteria
-        X = check_array(X,
-                        force_all_finite='allow-nan',
-                        estimator=self,
-                        copy=True)
+        X = check_array(X, force_all_finite="allow-nan", estimator=self, copy=True)
 
         # check to see if there are any null values, and if
         # so impute using the desired imputation approach
@@ -685,9 +696,11 @@ def fit(self, X, y=None):
             self.mean_ = np.mean(X, axis=0)
             cov_mtx = cov(X)
         else:
-            error_msg = ('If `is_cov_matrix=True`, then the rows and column in the data '
-                         'set must be equal, and must equal the number of variables '
-                         'in your model.')
+            error_msg = (
+                "If `is_cov_matrix=True`, then the rows and column in the data "
+                "set must be equal, and must equal the number of variables "
+                "in your model."
+            )
             assert X.shape[0] == X.shape[1] == self.model.n_variables, error_msg
             cov_mtx = X.copy()
 
@@ -701,51 +714,55 @@ def fit(self, X, y=None):
         factor_covs_init = np.full((self.model.n_lower_diag, 1), 0.05)
 
         # we merge all of the arrays into a single 1d vector
-        x0 = self._combine(loading_init,
-                           error_vars_init,
-                           factor_vars_init,
-                           factor_covs_init,
-                           self.model.n_factors,
-                           self.model.n_variables,
-                           self.model.n_lower_diag)
+        x0 = self._combine(
+            loading_init,
+            error_vars_init,
+            factor_vars_init,
+            factor_covs_init,
+            self.model.n_factors,
+            self.model.n_variables,
+            self.model.n_lower_diag,
+        )
 
         # if the bounds argument is None, then we initialized the
         # boundaries to (None, None) for everything except factor covariances;
         # at some point in the future, we may update this to place limits
         # on the loading matrix boundaries, too, but the case in R and SAS
         if self.bounds is not None:
-            error_msg = ('The length of `bounds` must equal the length of your '
-                         'input array `x0`: {} != {}.'.format(len(self.bounds), len(x0)))
+            error_msg = (
+                "The length of `bounds` must equal the length of your "
+                "input array `x0`: {} != {}.".format(len(self.bounds), len(x0))
+            )
             assert len(self.bounds) == len(x0), error_msg
 
         # fit the actual model using L-BFGS algorithm;
         # the constraints are set inside the objective function,
         # so that we can avoid using linear programming methods (e.g. SLSQP)
-        res = minimize(self._objective, x0,
-                       method='L-BFGS-B',
-                       options={'maxiter': self.max_iter, 'disp': self.disp},
-                       bounds=self.bounds,
-                       args=(cov_mtx, self.model.loadings))
+        res = minimize(
+            self._objective,
+            x0,
+            method="L-BFGS-B",
+            options={"maxiter": self.max_iter, "disp": self.disp},
+            bounds=self.bounds,
+            args=(cov_mtx, self.model.loadings),
+        )
 
         # if the optimizer failed to converge, print the message
         if not res.success:
-            warnings.warn('The optimization routine failed '
-                          'to converge: {}'.format(str(res.message)))
+            warnings.warn(
+                f"The optimization routine failed to converge: {str(res.message)}"
+            )
 
         # we split all the 1d array back into the set of original arrays
-        (loadings_res,
-         error_vars_res,
-         factor_vars_res,
-         factor_covs_res) = self._split(res.x,
-                                        self.model.n_factors,
-                                        self.model.n_variables,
-                                        self.model.n_lower_diag)
+        (loadings_res, error_vars_res, factor_vars_res, factor_covs_res) = self._split(
+            res.x, self.model.n_factors, self.model.n_variables, self.model.n_lower_diag
+        )
 
         # we combine the factor covariances and variances into
         # a single variance-covariance matrix to make things easier,
         # but also check to make see if anything was fixed
         factor_varcovs_res = merge_variance_covariance(factor_vars_res, factor_covs_res)
-        with np.errstate(all='ignore'):
+        with np.errstate(all="ignore"):
             factor_varcovs_res = covariance_to_correlation(factor_varcovs_res)
 
         self.loadings_ = loadings_res
@@ -756,21 +773,23 @@ def fit(self, X, y=None):
         self.log_likelihood_ = -res.fun
         self.aic_ = 2 * res.fun + 2 * (x0.shape[0] + self.model.n_variables)
         if self.n_obs is not None:
-            self.bic_ = 2 * res.fun + np.log(self.n_obs) * (x0.shape[0] + self.model.n_variables)
+            self.bic_ = 2 * res.fun + np.log(self.n_obs) * (
+                x0.shape[0] + self.model.n_variables
+            )
         return self
 
     def transform(self, X):
         """
-        Get the factor scores for new data set.
+        Get the factor scores for a new data set.
 
         Parameters
         ----------
-        X : array-like, shape (n_samples, n_features)
+        X : array-like, shape (``n_samples``, ``n_features``)
             The data to score using the fitted factor model.
 
         Returns
         -------
-        scores : numpy array, shape (n_samples, n_components)
+        scores : numpy array, shape (``n_samples``, ``n_components``)
             The latent variables of X.
 
         Examples
@@ -797,7 +816,6 @@ def transform(self, X):
         ----------
         https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6157408/
         """
-
         # check if the data is a data frame,
         # so we can convert it to an array
         if isinstance(X, pd.DataFrame):
@@ -805,24 +823,23 @@ def transform(self, X):
 
         # now check the array, and make sure it
         # meets all of our expected criteria
-        X = check_array(X,
-                        force_all_finite=True,
-                        estimator=self,
-                        copy=True)
+        X = check_array(X, force_all_finite=True, estimator=self, copy=True)
 
         # meets all of our expected criteria
-        check_is_fitted(self, ['loadings_', 'error_vars_'])
+        check_is_fitted(self, ["loadings_", "error_vars_"])
 
         # see if we saved the original mean and std
         if self.mean_ is None:
-            warnings.warn('Could not find original mean; using the mean '
-                          'from the current data set.')
+            warnings.warn(
+                "Could not find original mean; using the mean "
+                "from the current data set."
+            )
             mean = np.mean(X, axis=0)
         else:
             mean = self.mean_
 
         # get the scaled data
-        X_scale = (X - mean)
+        X_scale = X - mean
 
         # get the loadings and error variances
         loadings = self.loadings_.copy()
@@ -836,23 +853,23 @@ def transform(self, X):
 
         # calculate the weights, using Bartlett's formula
         # (lambda' error_covsˆ−1 lambda)ˆ−1 lambda' error_covsˆ−1 (X - muX)
-        weights = np.linalg.pinv(loadings.T.dot(error_covs_inv)
-                                           .dot(loadings)) \
-                           .dot(loadings.T).dot(error_covs_inv)
+        weights = (
+            np.linalg.pinv(loadings.T.dot(error_covs_inv).dot(loadings))
+            .dot(loadings.T)
+            .dot(error_covs_inv)
+        )
 
         scores = weights.dot(X_scale.T).T
         return scores
 
     def get_model_implied_cov(self):
         """
-        Get the model-implied covariance
-        matrix (sigma), if the model has been estimated.
+        Get the model-implied covariance matrix (sigma) for an estimated model.
 
         Returns
         -------
-        model_implied_cov : numpy array
-            The model-implied covariance
-            matrix.
+        model_implied_cov : np.ndarray
+            The model-implied covariance matrix.
 
         Examples
         --------
@@ -884,26 +901,25 @@ def get_model_implied_cov(self):
                 0.36730952, 0.1882011 , 1.28888233]])
         """
         # meets all of our expected criteria
-        check_is_fitted(self, ['loadings_', 'factor_varcovs_'])
+        check_is_fitted(self, ["loadings_", "factor_varcovs_"])
         error = np.diag(self.error_vars_.flatten())
         return self.loadings_.dot(self.factor_varcovs_).dot(self.loadings_.T) + error
 
     def _get_derivatives_implied_cov(self):
         """
-        Compute the derivatives for the implied covariance
-        matrix (sigma).
+        Compute the derivatives for the implied covariance matrix (sigma).
 
         Returns
         -------
-        loadings_dx : numpy array
+        loadings_dx : np.ndarray
             The derivative of the loadings matrix.
-        factor_covs_dx : numpy array
+        factor_covs_dx : np.ndarray
             The derivative of the factor covariance matrix.
-        error_covs_dx : numpy array
+        error_covs_dx : np.ndarray
             The derivative of the error covariance matrix.
         """
         # meets all of our expected criteria
-        check_is_fitted(self, 'loadings_')
+        check_is_fitted(self, "loadings_")
 
         loadings = self.loadings_.copy()
         factor_covs = self.factor_varcovs_.copy()
@@ -913,25 +929,31 @@ def _get_derivatives_implied_cov(self):
         sym_lower_fac_idx = get_symmetric_lower_idxs(self.model.n_factors, diag=False)
 
         factors_diag = np.eye(self.model.n_factors)
-        factors_diag_mult = factors_diag.dot(factor_covs).dot(factors_diag.T).dot(loadings.T)
+        factors_diag_mult = (
+            factors_diag.dot(factor_covs).dot(factors_diag.T).dot(loadings.T)
+        )
 
         # calculate the derivative of the loadings matrix, using the commutation matrix
-        loadings_dx = np.eye(self.model.n_variables**2) + commutation_matrix(self.model.n_variables,
-                                                                             self.model.n_variables)
-        loadings_dx = loadings_dx.dot(np.kron(factors_diag_mult, np.eye(self.model.n_variables)).T)
+        loadings_dx = np.eye(self.model.n_variables**2) + commutation_matrix(
+            self.model.n_variables, self.model.n_variables
+        )
+        loadings_dx = loadings_dx.dot(
+            np.kron(factors_diag_mult, np.eye(self.model.n_variables)).T
+        )
 
         # calculate the derivative of the factor_covs matrix
         factor_covs_dx = loadings.dot(factors_diag)
         factor_covs_dx = np.kron(factor_covs_dx, factor_covs_dx)
 
-        off_diag = (factor_covs_dx[:, sym_lower_fac_idx] +
-                    factor_covs_dx[:, sym_upper_fac_idx])
+        off_diag = (
+            factor_covs_dx[:, sym_lower_fac_idx] + factor_covs_dx[:, sym_upper_fac_idx]
+        )
 
         combine_indices = np.concatenate([sym_upper_fac_idx, sym_lower_fac_idx])
         combine_diag = np.concatenate([off_diag, off_diag], axis=1)
 
         factor_covs_dx[:, combine_indices] = combine_diag
-        factor_covs_dx = factor_covs_dx[:, :factor_covs.size]
+        factor_covs_dx = factor_covs_dx[:, : factor_covs.size]
 
         # calculate the derivative of the error_cov matrix,
         # which we assume will always be a diagonal matrix
@@ -943,33 +965,37 @@ def _get_derivatives_implied_cov(self):
         error_covs_dx = error_covs_dx[sym_lower_var_idx, :]
 
         # we also want to calculate the derivative for the intercepts
-        intercept_dx = np.zeros((loadings_dx.shape[0], self.model.n_variables), dtype=float)
-        return (loadings_dx[:, self.model.loadings_free].copy(),
-                factor_covs_dx[:, self.model.factor_covs_free].copy(),
-                error_covs_dx[:, self.model.error_vars_free].copy(),
-                intercept_dx)
+        intercept_dx = np.zeros(
+            (loadings_dx.shape[0], self.model.n_variables), dtype=float
+        )
+        return (
+            loadings_dx[:, self.model.loadings_free].copy(),
+            factor_covs_dx[:, self.model.factor_covs_free].copy(),
+            error_covs_dx[:, self.model.error_vars_free].copy(),
+            intercept_dx,
+        )
 
     def _get_derivatives_implied_mu(self):
         """
         Compute the "derivatives" for the implied means.
-        Note that the derivatives of the implied means
-        will not correspond to the actual mean values
-        of the original data set, because that data could be
-        a covariance matrix, rather than the full data set.
-        Thus, we assume the mean values are zero and the
-        data are normally distributed.
+
+        Note that the derivatives of the implied means will not correspond
+        to the actual mean values of the original data set, because that
+        data could be a covariance matrix, rather than the full data set.
+        Thus, we assume the mean values are zero and the data are normally
+        distributed.
 
         Returns
         -------
-        loadings_dx : numpy array
+        loadings_dx : np.ndarray
             The derivative of the loadings means.
-        factor_covs_dx : numpy array
+        factor_covs_dx : np.ndarray
             The derivative of the factor covariance means.
-        error_covs_dx : numpy array
+        error_covs_dx : np.ndarray
             The derivative of the error covariance means.
         """
         # meets all of our expected criteria
-        check_is_fitted(self, 'loadings_')
+        check_is_fitted(self, "loadings_")
 
         # initialize some matrices that we'll use below to
         # correct the shape of the mean loadings derivatives
@@ -979,34 +1005,37 @@ def _get_derivatives_implied_mu(self):
         # the mean derivatives will just be zeros for both
         # the error covariance and factor covariance matrices,
         # since we don't have actual mean values
-        error_covs_dx = np.zeros((self.model.n_variables, len(self.model.error_vars_free)))
-        factor_covs_dx = np.zeros((self.model.n_variables, len(self.model.factor_covs_free)))
+        error_covs_dx = np.zeros(
+            (self.model.n_variables, len(self.model.error_vars_free))
+        )
+        factor_covs_dx = np.zeros(
+            (self.model.n_variables, len(self.model.factor_covs_free))
+        )
 
         # again, the implied means are just going to be diagonal matrices
         # corresponding to the number of variables and factors
-        loadings_dx = np.kron(factors_diag.dot(factors_zero).T, np.eye(self.model.n_variables))
+        loadings_dx = np.kron(
+            factors_diag.dot(factors_zero).T, np.eye(self.model.n_variables)
+        )
         loadings_dx = loadings_dx[:, self.model.loadings_free].copy()
 
         # we also calculate the derivative for the intercept, which will be zeros again
         intercept_dx = np.zeros((loadings_dx.shape[0], self.model.n_variables))
-        intercept_dx[:self.model.n_variables,
-                     :self.model.n_variables] = np.eye(self.model.n_variables)
+        intercept_dx[: self.model.n_variables, : self.model.n_variables] = np.eye(
+            self.model.n_variables
+        )
 
-        return (loadings_dx,
-                factor_covs_dx,
-                error_covs_dx,
-                intercept_dx)
+        return (loadings_dx, factor_covs_dx, error_covs_dx, intercept_dx)
 
     def get_standard_errors(self):
         """
-        Get the standard errors from the implied
-        covariance matrix and implied means.
+        Get standard errors from the implied covariance matrix and implied means.
 
         Returns
         -------
-        loadings_se : numpy array
+        loadings_se : np.ndarray
             The standard errors for the factor loadings.
-        error_vars_se : numpy array
+        error_vars_se : np.ndarray
             The standard errors for the error variances.
 
         Examples
@@ -1033,17 +1062,21 @@ def get_standard_errors(self):
                 0.10176889, 0.06607524, 0.04742321, 0.05373646]))
         """
         # meets all of our expected criteria
-        check_is_fitted(self, ['loadings_', 'n_obs'])
-
-        (loadings_dx,
-         factor_covs_dx,
-         error_covs_dx,
-         intercept_dx) = self._get_derivatives_implied_cov()
-
-        (loadings_dx_mu,
-         factor_covs_dx_mu,
-         error_covs_dx_mu,
-         intercept_dx_mu) = self._get_derivatives_implied_mu()
+        check_is_fitted(self, ["loadings_", "n_obs"])
+
+        (
+            loadings_dx,
+            factor_covs_dx,
+            error_covs_dx,
+            intercept_dx,
+        ) = self._get_derivatives_implied_cov()
+
+        (
+            loadings_dx_mu,
+            factor_covs_dx_mu,
+            error_covs_dx_mu,
+            intercept_dx_mu,
+        ) = self._get_derivatives_implied_mu()
 
         # combine all of our derivatives; below we will  merge all of these
         # together in a single matrix, delta, to use the delta rule
@@ -1065,10 +1098,9 @@ def get_standard_errors(self):
         h1_information = block_diag(sigma_inv, h1_information)
 
         # we concatenate all derivatives into a delta matrix
-        delta = np.concatenate((loadings_dx,
-                                error_cov_dx,
-                                factor_covs_dx,
-                                intercept_dx), axis=1)
+        delta = np.concatenate(
+            (loadings_dx, error_cov_dx, factor_covs_dx, intercept_dx), axis=1
+        )
 
         # calculate the fisher information matrix
         information = delta.T.dot(h1_information).dot(delta)
@@ -1091,8 +1123,10 @@ def get_standard_errors(self):
         # format of the original loadings matrix
         loadings_se = np.zeros((self.model.n_factors * self.model.n_variables,))
         loadings_se[self.model.loadings_free] = se[:loadings_idx]
-        loadings_se = loadings_se.reshape((self.model.n_variables, self.model.n_factors), order='F')
+        loadings_se = loadings_se.reshape(
+            (self.model.n_variables, self.model.n_factors), order="F"
+        )
 
         # get the error variance standard errors
-        error_vars_se = se[loadings_idx: error_vars_idx]
+        error_vars_se = se[loadings_idx:error_vars_idx]
         return loadings_se, error_vars_se
diff --git a/factor_analyzer/factor_analyzer.py b/factor_analyzer/factor_analyzer.py
index 5e25332..c1a004f 100644
--- a/factor_analyzer/factor_analyzer.py
+++ b/factor_analyzer/factor_analyzer.py
@@ -1,10 +1,10 @@
 """
 Factor analysis using MINRES or ML, with optional rotation using Varimax or Promax.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
 :author: Nitin Madnani (nmadnani@ets.org)
 :organization: Educational Testing Service
-:date: 2021-10-18
+:date: 2022-09-05
 """
 
 import warnings
@@ -22,21 +22,20 @@
 from .rotator import OBLIQUE_ROTATIONS, POSSIBLE_ROTATIONS, Rotator
 from .utils import corr, impute_values, partial_correlations, smc
 
-POSSIBLE_SVDS = ['randomized', 'lapack']
+POSSIBLE_SVDS = ["randomized", "lapack"]
 
-POSSIBLE_IMPUTATIONS = ['mean', 'median', 'drop']
+POSSIBLE_IMPUTATIONS = ["mean", "median", "drop"]
 
-POSSIBLE_METHODS = ['ml', 'mle', 'uls', 'minres', 'principal']
+POSSIBLE_METHODS = ["ml", "mle", "uls", "minres", "principal"]
 
 
 def calculate_kmo(x):
     """
-    Calculate the Kaiser-Meyer-Olkin criterion
-    for items and overall. This statistic represents
-    the degree to which each observed variable is
-    predicted, without error, by the other variables
-    in the dataset. In general, a KMO < 0.6 is considered
-    inadequate.
+    Calculate the Kaiser-Meyer-Olkin criterion for items and overall.
+
+    This statistic represents the degree to which each observed variable is
+    predicted, without error, by the other variables in the dataset.
+    In general, a KMO < 0.6 is considered inadequate.
 
     Parameters
     ----------
@@ -45,12 +44,11 @@ def calculate_kmo(x):
 
     Returns
     -------
-    kmo_per_variable : numpy array
+    kmo_per_variable : np.ndarray
         The KMO score per item.
     kmo_total : float
-        The KMO score overall.
+        The overall KMO score.
     """
-
     # calculate the partial correlations
     partial_corr = partial_correlations(x)
 
@@ -79,8 +77,7 @@ def calculate_kmo(x):
 
 def calculate_bartlett_sphericity(x):
     """
-    Test the hypothesis that the correlation matrix
-    is equal to the identity matrix.identity
+    Compute the Bartlett sphericity test.
 
     H0: The matrix of population correlations is equal to I.
     H1: The matrix of population correlations is not equal to I.
@@ -95,7 +92,7 @@ def calculate_bartlett_sphericity(x):
     Parameters
     ----------
     x : array-like
-        The array from which to calculate sphericity.
+        The array for which to calculate sphericity.
 
     Returns
     -------
@@ -116,7 +113,9 @@ def calculate_bartlett_sphericity(x):
 
 class FactorAnalyzer(BaseEstimator, TransformerMixin):
     """
-    A FactorAnalyzer class, which -
+    The main exploratory factor analysis class.
+
+    This class:
         (1) Fits a factor analysis model using minres, maximum likelihood,
             or principal factor extraction and returns the loading matrix
         (2) Optionally performs a rotation, with method including:
@@ -135,12 +134,11 @@ class FactorAnalyzer(BaseEstimator, TransformerMixin):
         The number of factors to select.
         Defaults to 3.
     rotation : str, optional
-        The type of rotation to perform after
-        fitting the factor analysis model.
-        If set to None, no rotation will be performed,
-        nor will any associated Kaiser normalization.
+        The type of rotation to perform after fitting the factor analysis
+        model. If set to ``None``, no rotation will be performed, nor will
+        any associated Kaiser normalization.
 
-        Methods include:
+        Possible values include:
 
             (a) varimax (orthogonal rotation)
             (b) promax (oblique rotation)
@@ -153,64 +151,53 @@ class FactorAnalyzer(BaseEstimator, TransformerMixin):
         Defaults to 'promax'.
 
     method : {'minres', 'ml', 'principal'}, optional
-        The fitting method to use, either MINRES or
-        Maximum Likelihood.
+        The fitting method to use, either MINRES or Maximum Likelihood.
         Defaults to 'minres'.
     use_smc : bool, optional
-        Whether to use squared multiple correlation
-        as starting guesses for factor analysis.
-        Defaults to True.
+        Whether to use squared multiple correlation as starting guesses for
+        factor analysis.
+        Defaults to ``True``.
     bounds : tuple, optional
-        The lower and upper bounds on the variables
-        for "L-BFGS-B" optimization.
+        The lower and upper bounds on the variables for "L-BFGS-B" optimization.
         Defaults to (0.005, 1).
     impute : {'drop', 'mean', 'median'}, optional
-        If missing values are present in the data, either use
-        list-wise deletion ('drop') or impute the column median
-        ('median') or column mean ('mean').
+        How to handle missing values, if any, in the data: (a) use list-wise
+        deletion ('drop'), or (b) impute the column median ('median'), or
+        impute the column mean ('mean').
         Defaults to 'median'
-    use_corr_matrix : bool, optional
-        Set to true if the `data` is the correlation
-        matrix.
-        Defaults to False.
+    is_corr_matrix : bool, optional
+        Set to ``True if the ``data`` is the correlation matrix.
+        Defaults to `False`.
     svd_method : {‘lapack’, ‘randomized’}
-        The SVD method to use when ``method='principal'``.
-        If 'lapack', use standard SVD from ``scipy.linalg``.
-        If 'randomized', use faster ``randomized_svd``
-        function from scikit-learn. The latter should only
-        be used if the number of columns is greater than or
-        equal to the number of rows in in the dataset.
+        The SVD method to use when ``method`` is 'principal'. If 'lapack',
+        use standard SVD from ``scipy.linalg``. If 'randomized', use faster
+        ``randomized_svd`` function from scikit-learn. The latter should only
+        be used if the number of columns is greater than or equal to the
+        number of rows in in the dataset.
         Defaults to 'randomized'
     rotation_kwargs, optional
-        Additional key word arguments
-        are passed to the rotation method.
+        Dictionary containing keyword arguments for the rotation method.
 
     Attributes
     ----------
     loadings : numpy array
         The factor loadings matrix.
-        Default to None, if `fit()` has not
-        been called.
+        None, if ``fit()``` has not been called.
     corr : numpy array
         The original correlation matrix.
-        Default to None, if `fit()` has not
-        been called.
+        ``None``, if ``fit()`` has not been called.
     rotation_matrix : numpy array
-        The rotation matrix, if a rotation
-        has been performed.
+        The rotation matrix, if a rotation has been performed.
     structure :numpy array or None
-        The structure loading matrix.
-        This only exists if the rotation
-        is promax.
+        The structure loading matrix. This only exists if ``rotation``
+        is 'promax'.
     psi : numpy array or None
-        The factor correlations
-        matrix. This only exists
-        if the rotation is oblique.
+        The factor correlations matrix. This only exists if ``rotation``
+        is 'oblique'.
 
     Notes
     -----
-    This code was partly derived from the excellent R package
-    `psych`.
+    This code was partly derived from the excellent R package `psych`.
 
     References
     ----------
@@ -242,17 +229,19 @@ class FactorAnalyzer(BaseEstimator, TransformerMixin):
            0.66208428, 0.61911036, 0.73194557, 0.64929612, 0.71149718])
     """
 
-    def __init__(self,
-                 n_factors=3,
-                 rotation='promax',
-                 method='minres',
-                 use_smc=True,
-                 is_corr_matrix=False,
-                 bounds=(0.005, 1),
-                 impute='median',
-                 svd_method='randomized',
-                 rotation_kwargs=None):
-
+    def __init__(
+        self,
+        n_factors=3,
+        rotation="promax",
+        method="minres",
+        use_smc=True,
+        is_corr_matrix=False,
+        bounds=(0.005, 1),
+        impute="median",
+        svd_method="randomized",
+        rotation_kwargs=None,
+    ):
+        """Initialize the factor analyzer."""
         self.n_factors = n_factors
         self.rotation = rotation
         self.method = method
@@ -278,35 +267,58 @@ def __init__(self,
     def _arg_checker(self):
         """
         Check the input parameters to make sure they're properly formattted.
-        We need to do this to ensure that the FactorAnalyzer class can be properly
-        cloned when used with grid search CV, for example.
+
+        We need to do this to ensure that the FactorAnalyzer class can be
+        properly cloned when used with grid search CV, for example.
         """
-        self.rotation = self.rotation.lower() if isinstance(self.rotation, str) else self.rotation
+        self.rotation = (
+            self.rotation.lower() if isinstance(self.rotation, str) else self.rotation
+        )
         if self.rotation not in POSSIBLE_ROTATIONS + [None]:
-            raise ValueError(f"The rotation must be one of the following: {POSSIBLE_ROTATIONS + [None]}")
+            raise ValueError(
+                f"The rotation must be one of the following: {POSSIBLE_ROTATIONS + [None]}"
+            )
 
-        self.method = self.method.lower() if isinstance(self.method, str) else self.method
+        self.method = (
+            self.method.lower() if isinstance(self.method, str) else self.method
+        )
         if self.method not in POSSIBLE_METHODS:
-            raise ValueError(f"The method must be one of the following: {POSSIBLE_METHODS}")
+            raise ValueError(
+                f"The method must be one of the following: {POSSIBLE_METHODS}"
+            )
 
-        self.impute = self.impute.lower() if isinstance(self.impute, str) else self.impute
+        self.impute = (
+            self.impute.lower() if isinstance(self.impute, str) else self.impute
+        )
         if self.impute not in POSSIBLE_IMPUTATIONS:
-            raise ValueError(f"The imputation must be one of the following: {POSSIBLE_IMPUTATIONS}")
-
-        self.svd_method = self.svd_method.lower() if isinstance(self.svd_method, str) else self.svd_method
+            raise ValueError(
+                f"The imputation must be one of the following: {POSSIBLE_IMPUTATIONS}"
+            )
+
+        self.svd_method = (
+            self.svd_method.lower()
+            if isinstance(self.svd_method, str)
+            else self.svd_method
+        )
         if self.svd_method not in POSSIBLE_SVDS:
-            raise ValueError(f"The SVD method must be one of the following: {POSSIBLE_SVDS}")
+            raise ValueError(
+                f"The SVD method must be one of the following: {POSSIBLE_SVDS}"
+            )
 
-        if self.method == 'principal' and self.is_corr_matrix:
-            raise ValueError('The principal method is only implemented using '
-                             'the full data set, not the correlation matrix.')
+        if self.method == "principal" and self.is_corr_matrix:
+            raise ValueError(
+                "The principal method is only implemented using "
+                "the full data set, not the correlation matrix."
+            )
 
-        self.rotation_kwargs = {} if self.rotation_kwargs is None else self.rotation_kwargs
+        self.rotation_kwargs = (
+            {} if self.rotation_kwargs is None else self.rotation_kwargs
+        )
 
     @staticmethod
-    def _fit_uls_objective(psi, corr_mtx, n_factors):
+    def _fit_uls_objective(psi, corr_mtx, n_factors):  # noqa: D401
         """
-        The objective function passed to `minimize()` for ULS.
+        The objective function passed for unweighted least-squares (ULS).
 
         Parameters
         ----------
@@ -320,8 +332,8 @@ def _fit_uls_objective(psi, corr_mtx, n_factors):
         Returns
         -------
         error : float
-            The scalar error calculated from the residuals
-            of the loading matrix.
+            The scalar error calculated from the residuals of the loading
+            matrix.
         """
         np.fill_diagonal(corr_mtx, 1 - psi)
 
@@ -352,15 +364,14 @@ def _fit_uls_objective(psi, corr_mtx, n_factors):
         # lower triangle of the residual matrix; this could be
         # implemented here using `np.tril()` when this change is
         # merged into the stable version of `psych`.
-        residual = (corr_mtx - model)**2
+        residual = (corr_mtx - model) ** 2
         error = sp.sum(residual)
         return error
 
     @staticmethod
     def _normalize_uls(solution, corr_mtx, n_factors):
         """
-        Weighted least squares normalization for loadings
-        estimated using MINRES.
+        Weighted least squares normalization for loadings using MINRES.
 
         Parameters
         ----------
@@ -373,7 +384,7 @@ def _normalize_uls(solution, corr_mtx, n_factors):
 
         Returns
         -------
-        loadings : numpy array
+        loadings : np.ndarray
             The factor loading matrix
         """
         np.fill_diagonal(corr_mtx, 1 - solution)
@@ -391,9 +402,9 @@ def _normalize_uls(solution, corr_mtx, n_factors):
         return loadings
 
     @staticmethod
-    def _fit_ml_objective(psi, corr_mtx, n_factors):
+    def _fit_ml_objective(psi, corr_mtx, n_factors):  # noqa: D401
         """
-        The objective function passed to `minimize()` for ML.
+        The objective function for maximum likelihood.
 
         Parameters
         ----------
@@ -412,9 +423,9 @@ def _fit_ml_objective(psi, corr_mtx, n_factors):
 
         Note
         ----
-        The ML objective is based on the `factanal()` function
-        from R's `stats` package. It may generate results different
-        from the `fa()` function in `psych`.
+        The ML objective is based on the `factanal()` function from ``stats``
+        package in R. It may generate different results from the ``fa()``
+        function in ``psych``.
 
         References
         ----------
@@ -428,14 +439,13 @@ def _fit_ml_objective(psi, corr_mtx, n_factors):
         values = values[::-1][n_factors:]
 
         # calculate the error
-        error = -(np.sum(np.log(values) - values) -
-                  n_factors + corr_mtx.shape[0])
+        error = -(np.sum(np.log(values) - values) - n_factors + corr_mtx.shape[0])
         return error
 
     @staticmethod
     def _normalize_ml(solution, corr_mtx, n_factors):
         """
-        Normalization for loadings estimated using ML.
+        Normalize loadings estimated using maximum likelihood.
 
         Parameters
         ----------
@@ -448,7 +458,7 @@ def _normalize_ml(solution, corr_mtx, n_factors):
 
         Returns
         -------
-        loadings : numpy array
+        loadings : np.ndarray
             The factor loading matrix
         """
         sc = np.diag(1 / np.sqrt(solution))
@@ -470,8 +480,7 @@ def _normalize_ml(solution, corr_mtx, n_factors):
 
     def _fit_principal(self, X):
         """
-        Fit the factor analysis model using a principal
-        factor analysis solution.
+        Fit factor analysis model using principal factor analysis.
 
         Parameters
         ----------
@@ -480,7 +489,7 @@ def _fit_principal(self, X):
 
         Returns
         -------
-        loadings : numpy array
+        loadings : np.ndarray
             The factor loadings matrix.
         """
         # standardize the data
@@ -491,16 +500,18 @@ def _fit_principal(self, X):
         # warn the user that the number of factors will be constrained
         nrows, ncols = X.shape
         if nrows < ncols and self.n_factors >= nrows:
-            warnings.warn('The number of factors will be '
-                          'constrained to min(n_samples, n_features)'
-                          '={}.'.format(min(nrows, ncols)))
+            warnings.warn(
+                "The number of factors will be "
+                "constrained to min(n_samples, n_features)"
+                "={}.".format(min(nrows, ncols))
+            )
 
         # perform the randomized singular value decomposition
-        if self.svd_method == 'randomized':
-            U, S, V = randomized_svd(X, self.n_factors, random_state=0)
+        if self.svd_method == "randomized":
+            _, _, V = randomized_svd(X, self.n_factors, random_state=0)
         # otherwise, perform the full SVD
         else:
-            U, S, V = np.linalg.svd(X, full_matrices=False)
+            _, _, V = np.linalg.svd(X, full_matrices=False)
 
         corr_mtx = np.dot(X, V.T)
         loadings = np.array([[pearsonr(x, c)[0] for c in corr_mtx.T] for x in X.T])
@@ -508,8 +519,7 @@ def _fit_principal(self, X):
 
     def _fit_factor_analysis(self, corr_mtx):
         """
-        Fit the factor analysis model using either
-        minres or ml solutions.
+        Fit factor analysis model using either MINRES or maximum likelihood.
 
         Parameters
         ----------
@@ -518,16 +528,14 @@ def _fit_factor_analysis(self, corr_mtx):
 
         Returns
         -------
-        loadings : numpy array
-            The factor loadings matrix.
+        loadings : np.ndarray
 
         Raises
         ------
         ValueError
-            If any of the correlations are null, most likely due
-            to having zero standard deviation.
+            If any of the correlations are null, most likely due to having
+            zero standard deviation.
         """
-
         # if `use_smc` is True, get get squared multiple correlations
         # and use these as initial guesses for optimizer
         if self.use_smc:
@@ -546,35 +554,37 @@ def _fit_factor_analysis(self, corr_mtx):
 
         # minimize the appropriate objective function
         # and the L-BFGS-B algorithm
-        if self.method == 'ml' or self.method == 'mle':
+        if self.method == "ml" or self.method == "mle":
             objective = self._fit_ml_objective
         else:
             objective = self._fit_uls_objective
 
         # use scipy to perform the actual minimization
-        res = minimize(objective,
-                       start,
-                       method='L-BFGS-B',
-                       bounds=bounds,
-                       options={'maxiter': 1000},
-                       args=(corr_mtx, self.n_factors))
+        res = minimize(
+            objective,
+            start,
+            method="L-BFGS-B",
+            bounds=bounds,
+            options={"maxiter": 1000},
+            args=(corr_mtx, self.n_factors),
+        )
 
         if not res.success:
-            warnings.warn('Failed to converge: {}'.format(res.message))
+            warnings.warn(f"Failed to converge: {res.message}")
 
         # transform the final loading matrix (using wls for MINRES,
         # and ml normalization for ML), and convert to DataFrame
-        if self.method == 'ml' or self.method == 'mle':
-            loadings = self._normalize_ml(res.x, corr_mtx, self. n_factors)
+        if self.method == "ml" or self.method == "mle":
+            loadings = self._normalize_ml(res.x, corr_mtx, self.n_factors)
         else:
             loadings = self._normalize_uls(res.x, corr_mtx, self.n_factors)
         return loadings
 
     def fit(self, X, y=None):
         """
-        Fit the factor analysis model using either
-        minres, ml, or principal solutions. By default, use SMC
-        as starting guesses.
+        Fit factor analysis model using either MINRES, ML, or principal factor analysis.
+
+        By default, use SMC as starting guesses.
 
         Parameters
         ----------
@@ -604,7 +614,6 @@ def fit(self, X, y=None):
                [ 0.76041819, -0.23768727, -0.1206858 ],
                [ 0.81533404, -0.12494695,  0.17639683]])
         """
-
         # check the input arguments
         self._arg_checker()
 
@@ -617,10 +626,7 @@ def fit(self, X, y=None):
 
         # now check the array, and make sure it
         # meets all of our expected criteria
-        X = check_array(X,
-                        force_all_finite='allow-nan',
-                        estimator=self,
-                        copy=True)
+        X = check_array(X, force_all_finite="allow-nan", estimator=self, copy=True)
 
         # check to see if there are any null values, and if
         # so impute using the desired imputation approach
@@ -639,7 +645,7 @@ def fit(self, X, y=None):
         self.corr_ = corr_mtx.copy()
 
         # fit factor analysis model
-        if self.method == 'principal':
+        if self.method == "principal":
             loadings = self._fit_principal(X)
         else:
             loadings = self._fit_factor_analysis(corr_mtx)
@@ -653,19 +659,23 @@ def fit(self, X, y=None):
         # whether to rotate the loadings matrix
         if self.rotation is not None:
             if loadings.shape[1] <= 1:
-                warnings.warn('No rotation will be performed when '
-                              'the number of factors equals 1.')
+                warnings.warn(
+                    "No rotation will be performed when "
+                    "the number of factors equals 1."
+                )
             else:
-                if 'method' in self.rotation_kwargs:
-                    warnings.warn('You cannot pass a rotation method to '
-                                  '`rotation_kwargs`. This will be ignored.')
-                    self.rotation_kwargs.pop('method')
+                if "method" in self.rotation_kwargs:
+                    warnings.warn(
+                        "You cannot pass a rotation method to "
+                        "`rotation_kwargs`. This will be ignored."
+                    )
+                    self.rotation_kwargs.pop("method")
                 rotator = Rotator(method=self.rotation, **self.rotation_kwargs)
                 loadings = rotator.fit_transform(loadings)
                 rotation_mtx = rotator.rotation_
                 phi = rotator.phi_
                 # update the rotation matrix for everything, except promax
-                if self.rotation != 'promax':
+                if self.rotation != "promax":
                     rotation_mtx = np.linalg.inv(rotation_mtx).T
 
         if self.n_factors > 1:
@@ -679,11 +689,15 @@ def fit(self, X, y=None):
                 # update phi, if it exists -- that is, if the rotation is oblique
                 # create the structure matrix for any oblique rotation
                 phi = np.dot(np.dot(np.diag(signs), phi), np.diag(signs))
-                structure = np.dot(loadings, phi) if self.rotation in OBLIQUE_ROTATIONS else None
+                structure = (
+                    np.dot(loadings, phi)
+                    if self.rotation in OBLIQUE_ROTATIONS
+                    else None
+                )
 
         # resort the factors according to their variance,
         # unless the method is principal
-        if self.method != 'principal':
+        if self.method != "principal":
             variance = self._get_factor_variance(loadings)[0]
             new_order = list(reversed(np.argsort(variance)))
             loadings = loadings[:, new_order].copy()
@@ -701,16 +715,16 @@ def fit(self, X, y=None):
 
     def transform(self, X):
         """
-        Get the factor scores for new data set.
+        Get factor scores for a new data set.
 
         Parameters
         ----------
-        X : array-like, shape (n_samples, n_features)
+        X : array-like, shape (``n_samples``, ``n_features``)
             The data to score using the fitted factor model.
 
         Returns
         -------
-        X_new : numpy array, shape (n_samples, n_components)
+        X_new : np.ndarray, shape (``n_samples``, ``n_components``)
             The latent variables of X.
 
         Examples
@@ -732,7 +746,6 @@ def transform(self, X):
                [ 1.86904519, -0.3532394 , -0.68170573],
                [ 0.86133386,  0.18280695, -0.79170903]])
         """
-
         # check if the data is a data frame,
         # so we can convert it to an array
         if isinstance(X, pd.DataFrame):
@@ -742,18 +755,17 @@ def transform(self, X):
 
         # now check the array, and make sure it
         # meets all of our expected criteria
-        X = check_array(X,
-                        force_all_finite=True,
-                        estimator=self,
-                        copy=True)
+        X = check_array(X, force_all_finite=True, estimator=self, copy=True)
 
         # meets all of our expected criteria
-        check_is_fitted(self, 'loadings_')
+        check_is_fitted(self, "loadings_")
 
         # see if we saved the original mean and std
         if self.mean_ is None or self.std_ is None:
-            warnings.warn('Could not find original mean and standard deviation; using'
-                          'the mean and standard deviation from the current data set.')
+            warnings.warn(
+                "Could not find original mean and standard deviation; using"
+                "the mean and standard deviation from the current data set."
+            )
             mean = np.mean(X, axis=0)
             std = np.std(X, axis=0)
         else:
@@ -773,8 +785,10 @@ def transform(self, X):
         try:
             self.weights_ = np.linalg.solve(self.corr_, structure)
         except Exception as error:
-            warnings.warn('Unable to calculate the factor score weights; '
-                          'factor loadings used instead: {}'.format(error))
+            warnings.warn(
+                "Unable to calculate the factor score weights; "
+                "factor loadings used instead: {}".format(error)
+            )
             self.weights_ = self.loadings_
 
         scores = np.dot(X_scale, self.weights_)
@@ -782,15 +796,14 @@ def transform(self, X):
 
     def get_eigenvalues(self):
         """
-        Calculate the eigenvalues, given the
-        factor correlation matrix.
+        Calculate the eigenvalues, given the factor correlation matrix.
 
         Returns
         -------
-        original_eigen_values : numpy array
-            The original eigen values
-        common_factor_eigen_values : numpy array
-            The common factor eigen values
+        original_eigen_values : np.ndarray
+            The original eigenvalues
+        common_factor_eigen_values : np.ndarray
+            The common factor eigenvalues
 
         Examples
         --------
@@ -809,7 +822,7 @@ def get_eigenvalues(self):
                 0.01029184, -0.0074    , -0.03694834, -0.05959057, -0.07428059]))
         """
         # meets all of our expected criteria
-        check_is_fitted(self, ['loadings_', 'corr_'])
+        check_is_fitted(self, ["loadings_", "corr_"])
         corr_mtx = self.corr_.copy()
 
         e_values, _ = np.linalg.eigh(corr_mtx)
@@ -825,14 +838,12 @@ def get_eigenvalues(self):
 
     def get_communalities(self):
         """
-        Calculate the communalities, given the
-        factor loading matrix.
+        Calculate the communalities, given the factor loading matrix.
 
         Returns
         -------
-        communalities : numpy array
-            The communalities from the factor loading
-            matrix.
+        communalities : np.ndarray
+            The communalities from the factor loading matrix.
 
         Examples
         --------
@@ -849,21 +860,19 @@ def get_communalities(self):
                0.66208428, 0.61911036, 0.73194557, 0.64929612, 0.71149718])
         """
         # meets all of our expected criteria
-        check_is_fitted(self, 'loadings_')
+        check_is_fitted(self, "loadings_")
         loadings = self.loadings_.copy()
-        communalities = (loadings ** 2).sum(axis=1)
+        communalities = (loadings**2).sum(axis=1)
         return communalities
 
     def get_uniquenesses(self):
         """
-        Calculate the uniquenesses, given the
-        factor loading matrix.
+        Calculate the uniquenesses, given the factor loading matrix.
 
         Returns
         -------
-        uniquenesses : numpy array
-            The uniquenesses from the factor loading
-            matrix.
+        uniquenesses : np.ndarray
+            The uniquenesses from the factor loading matrix.
 
         Examples
         --------
@@ -880,38 +889,38 @@ def get_uniquenesses(self):
                0.33791572, 0.38088964, 0.26805443, 0.35070388, 0.28850282])
         """
         # meets all of our expected criteria
-        check_is_fitted(self, 'loadings_')
+        check_is_fitted(self, "loadings_")
         communalities = self.get_communalities()
         communalities = communalities.copy()
-        uniqueness = (1 - communalities)
+        uniqueness = 1 - communalities
         return uniqueness
 
     @staticmethod
     def _get_factor_variance(loadings):
         """
-        A helper method to get the factor variances,
-        because sometimes we need them even before the
-        model is fitted.
+        Get the factor variances.
+
+        This is a private helper method to get the factor variances,
+        because sometimes we need them even before the model is fitted.
 
         Parameters
         ----------
         loadings : array-like
-            The factor loading matrix,
-            in whatever state.
+            The factor loading matrix, in whatever state.
 
         Returns
         -------
-        variance : numpy array
+        variance : np.ndarray
             The factor variances.
-        proportional_variance : numpy array
+        proportional_variance : np.ndarray
             The proportional factor variances.
-        cumulative_variances : numpy array
+        cumulative_variances : np.ndarray
             The cumulative factor variances.
         """
         n_rows = loadings.shape[0]
 
         # calculate variance
-        loadings = loadings ** 2
+        loadings = loadings**2
         variance = np.sum(loadings, axis=0)
 
         # calculate proportional variance
@@ -920,23 +929,22 @@ def _get_factor_variance(loadings):
         # calculate cumulative variance
         cumulative_variance = np.cumsum(proportional_variance, axis=0)
 
-        return (variance,
-                proportional_variance,
-                cumulative_variance)
+        return (variance, proportional_variance, cumulative_variance)
 
     def get_factor_variance(self):
         """
-        Calculate the factor variance information,
-        including variance, proportional variance
-        and cumulative variance for each factor
+        Calculate factor variance information.
+
+        The factor variance information including the variance,
+        proportional variance, and cumulative variance for each factor.
 
         Returns
         -------
-        variance : numpy array
+        variance : np.ndarray
             The factor variances.
-        proportional_variance : numpy array
+        proportional_variance : np.ndarray
             The proportional factor variances.
-        cumulative_variances : numpy array
+        cumulative_variances : np.ndarray
             The cumulative factor variances.
 
         Examples
@@ -958,6 +966,6 @@ def get_factor_variance(self):
          array([0.35101885, 0.47938987, 0.55312938]))
         """
         # meets all of our expected criteria
-        check_is_fitted(self, 'loadings_')
+        check_is_fitted(self, "loadings_")
         loadings = self.loadings_.copy()
         return self._get_factor_variance(loadings)
diff --git a/factor_analyzer/rotator.py b/factor_analyzer/rotator.py
index 02b3bad..47dfc4a 100644
--- a/factor_analyzer/rotator.py
+++ b/factor_analyzer/rotator.py
@@ -1,28 +1,29 @@
 """
 Class to perform various rotations of factor loading matrices.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
 :author: Nitin Madnani (nmadnani@ets.org)
 :organization: Educational Testing Service
-:date: 2021-10-18
+:date: 2022-09-05
 """
 
 import numpy as np
 import scipy as sp
 from sklearn.base import BaseEstimator
 
-ORTHOGONAL_ROTATIONS = ['varimax', 'oblimax', 'quartimax', 'equamax', 'geomin_ort']
+ORTHOGONAL_ROTATIONS = ["varimax", "oblimax", "quartimax", "equamax", "geomin_ort"]
 
-OBLIQUE_ROTATIONS = ['promax', 'oblimin', 'quartimin', 'geomin_obl']
+OBLIQUE_ROTATIONS = ["promax", "oblimin", "quartimin", "geomin_obl"]
 
 POSSIBLE_ROTATIONS = ORTHOGONAL_ROTATIONS + OBLIQUE_ROTATIONS
 
 
 class Rotator(BaseEstimator):
     """
-    The Rotator class takes an (unrotated)
-    factor loading matrix and performs one
-    of several rotations.
+    Perform rotations on an unrotated factor loading matrix.
+
+    The Rotator class takes an (unrotated) factor loading matrix and
+    performs one of several rotations.
 
     Parameters
     ----------
@@ -37,54 +38,52 @@ class Rotator(BaseEstimator):
             (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
-        and de-normalization prior to and following
-        rotation. Used for varimax and promax rotations.
-        If None, default for promax is False, and default
-        for varimax is True.
-        Defaults to None.
+        Whether to perform Kaiser normalization and de-normalization prior
+        to and following rotation. Used for 'varimax' and 'promax' rotations.
+        If ``None``, default for 'promax' is ``False``, and default for
+        'varimax' is ``True``.
+        Defaults to ``None``.
     power : int, optional
-        The power to which to raise the promax loadings
-        (minus 1). Numbers should generally range form 2 to 4.
+        The exponent to which to raise the promax loadings (minus 1).
+        Numbers should generally range from 2 to 4.
         Defaults to 4.
-    kappa : int, optional
-        The kappa value for the equamax objective.
-        Ignored if the method is not 'equamax'.
+    kappa : float, optional
+        The kappa value for the 'equamax' objective. Ignored if the method
+        is not 'equamax'.
         Defaults to 0.
     gamma : int, optional
-        The gamma level for the oblimin objective.
-        Ignored if the method is not 'oblimin'.
+        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_*'.
+        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.
-        Defaults to `1000`.
+        The maximum number of iterations. Used for 'varimax' and 'oblique'
+        rotations.
+        Defaults to 1000.
     tol : float, optional
-        The convergence threshold.
-        Used for varimax and oblique rotations.
-        Defaults to `1e-5`.
+        The convergence threshold. Used for 'varimax' and 'oblique' rotations.
+        Defaults to 1e-5.
 
     Attributes
     ----------
-    loadings_ : numpy array, shape (n_features, n_factors)
-        The loadings matrix
-    rotation_ : numpy array, shape (n_factors, n_factors)
-        The rotation matrix
-    psi_ : numpy array or None
-        The factor correlations
-        matrix. This only exists
-        if the rotation is oblique.
+    loadings_ : np.ndarray, shape (``n_features``, ``n_factors``)
+        The loadings matrix.
+    rotation_ : np.ndarray, shape (``n_factors``, ``n_factors``)
+        The rotation matrix.
+    psi_ : np.ndarray or None
+        The factor correlations matrix. This only exists if ``method`` is
+        'oblique'.
 
     Notes
     -----
-    Most of the rotations in this class
-    are ported from R's `GPARotation` package.
+    Most of the rotations in this class are ported from R's ``GPARotation``
+    package.
 
     References
     ----------
@@ -111,16 +110,18 @@ class Rotator(BaseEstimator):
            [ 0.81372945,  0.20915845,  0.07479506]])
     """
 
-    def __init__(self,
-                 method='varimax',
-                 normalize=True,
-                 power=4,
-                 kappa=0,
-                 gamma=0,
-                 delta=0.01,
-                 max_iter=500,
-                 tol=1e-5):
-
+    def __init__(
+        self,
+        method="varimax",
+        normalize=True,
+        power=4,
+        kappa=0,
+        gamma=0,
+        delta=0.01,
+        max_iter=500,
+        tol=1e-5,
+    ):
+        """Initialize the rotator class."""
         self.method = method
         self.normalize = normalize
         self.power = power
@@ -134,7 +135,7 @@ def __init__(self,
         self.rotation_ = None
         self.phi_ = None
 
-    def _oblimax_obj(self, loadings):
+    def _oblimax_obj(self, loadings):  # noqa: D401
         """
         The Oblimax function objective.
 
@@ -146,16 +147,16 @@ def _oblimax_obj(self, loadings):
         Returns
         -------
         gradient_dict : dict
-            A dictionary with
-            - grad : np.array
-                The gradient.
-            - criterion : float
-                The value of the criterion for the objective.
+            A dictionary containing the following keys:
+            (1) ``grad`` : np.ndarray, containing the gradients.
+            (2) ``criterion`` : float, containing the criterion for the objective.
         """
-        gradient = -(4 * loadings**3 / (np.sum(loadings**4)) - 4 * loadings /
-                     (np.sum(loadings**2)))
-        criterion = (np.log(np.sum(loadings**4)) - 2 * np.log(np.sum(loadings**2)))
-        return {'grad': gradient, 'criterion': criterion}
+        gradient = -(
+            4 * loadings**3 / (np.sum(loadings**4))
+            - 4 * loadings / (np.sum(loadings**2))
+        )
+        criterion = np.log(np.sum(loadings**4)) - 2 * np.log(np.sum(loadings**2))
+        return {"grad": gradient, "criterion": criterion}
 
     def _quartimax_obj(self, loadings):
         """
@@ -164,22 +165,20 @@ def _quartimax_obj(self, loadings):
         Parameters
         ----------
         loadings : array-like
-            The loading matrix
+            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.
+            A dictionary containing the following keys:
+            (1) ``grad`` : np.ndarray, containing the gradients.
+            (2) ``criterion`` : float, containing the criterion for the objective.
         """
-        gradient = -loadings**3
+        gradient = -(loadings**3)
         criterion = -np.sum(np.diag(np.dot((loadings**2).T, loadings**2))) / 4
-        return {'grad': gradient, 'criterion': criterion}
+        return {"grad": gradient, "criterion": criterion}
 
-    def _oblimin_obj(self, loadings):
+    def _oblimin_obj(self, loadings):  # noqa: D401
         """
         The Oblimin function objective.
 
@@ -191,60 +190,54 @@ def _oblimin_obj(self, loadings):
         Returns
         -------
         gradient_dict : dict
-            A dictionary with
-            - grad : np.array
-                The gradient.
-            - criterion : float
-                The value of the criterion for the objective.
+            A dictionary containing the following keys:
+            (1) ``grad`` : np.ndarray, containing the gradients.
+            (2) ``criterion`` : float, containing the criterion for the objective.
         """
         X = np.dot(loadings**2, np.eye(loadings.shape[1]) != 1)
-        if (self.gamma != 0):
+        if self.gamma != 0:
             p = loadings.shape[0]
             X = np.diag(np.full(1, p)) - np.dot(np.zeros((p, p)), X)
         gradient = loadings * X
         criterion = np.sum(loadings**2 * X) / 4
-        return {'grad': gradient, 'criterion': criterion}
+        return {"grad": gradient, "criterion": criterion}
 
-    def _quartimin_obj(self, loadings):
+    def _quartimin_obj(self, loadings):  # noqa: D401
         """
-        Quartimin function objective.
+        The Quartimin function objective.
 
         Parameters
         ----------
         loadings : array-like
-            The loading matrix
+            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.
+            A dictionary containing the following keys:
+            (1) ``grad`` : np.ndarray, containing the gradients.
+            (2) ``criterion`` : float, containing the criterion for the objective.
         """
         X = np.dot(loadings**2, np.eye(loadings.shape[1]) != 1)
         gradient = loadings * X
         criterion = np.sum(loadings**2 * X) / 4
-        return {'grad': gradient, 'criterion': criterion}
+        return {"grad": gradient, "criterion": criterion}
 
-    def _equamax_obj(self, loadings):
+    def _equamax_obj(self, loadings):  # noqa: D401
         """
-        Equamax function objective.
+        The Equamax function objective.
 
         Parameters
         ----------
         loadings : array-like
-            The loading matrix
+            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.
+            A dictionary containing the following keys:
+            (1) ``grad`` : np.ndarray, containing the gradients.
+            (2) ``criterion`` : float, containing the criterion for the objective.
         """
         p, k = loadings.shape
 
@@ -252,34 +245,39 @@ def _equamax_obj(self, loadings):
         M = np.ones(p) - np.eye(p)
 
         loadings_squared = loadings**2
-        f1 = (1 - self.kappa) * np.sum(np.diag(np.dot(loadings_squared.T,
-                                                      np.dot(loadings_squared, N)))) / 4
-        f2 = self.kappa * np.sum(np.diag(np.dot(loadings_squared.T,
-                                                np.dot(M, loadings_squared)))) / 4
-
-        gradient = ((1 - self.kappa) * loadings * np.dot(loadings_squared, N) +
-                    self.kappa * loadings * np.dot(M, loadings_squared))
+        f1 = (
+            (1 - self.kappa)
+            * np.sum(np.diag(np.dot(loadings_squared.T, np.dot(loadings_squared, N))))
+            / 4
+        )
+        f2 = (
+            self.kappa
+            * np.sum(np.diag(np.dot(loadings_squared.T, np.dot(M, loadings_squared))))
+            / 4
+        )
+
+        gradient = (1 - self.kappa) * loadings * np.dot(
+            loadings_squared, N
+        ) + self.kappa * loadings * np.dot(M, loadings_squared)
 
         criterion = f1 + f2
-        return {'grad': gradient, 'criterion': criterion}
+        return {"grad": gradient, "criterion": criterion}
 
-    def _geomin_obj(self, loadings):
+    def _geomin_obj(self, loadings):  # noqa: D401
         """
-        Geomin function objective
+        The Geomin function objective.
 
         Parameters
         ----------
         loadings : array-like
-            The loading matrix
+            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.
+            A dictionary containing the following keys:
+            (1) ``grad`` : np.ndarray, containing the gradients.
+            (2) ``criterion`` : float, containing the criterion for the objective.
         """
         p, k = loadings.shape
 
@@ -290,14 +288,14 @@ def _geomin_obj(self, loadings):
 
         gradient = (2 / k) * (loadings / loadings2) * rep
         criterion = np.sum(pro)
-        return {'grad': gradient, 'criterion': criterion}
+        return {"grad": gradient, "criterion": criterion}
 
     def _oblique(self, loadings, method):
         """
-        A generic function for performing
-        all oblique rotations, except for
-        promax, which is implemented
-        separately.
+        Perform oblique rotations, except 'promax'.
+
+        A generic function for performing all oblique rotations, except for
+        promax, which is implemented separately.
 
         Parameters
         ----------
@@ -308,25 +306,23 @@ def _oblique(self, loadings, method):
 
         Returns
         -------
-        loadings : numpy array, shape (n_features, n_factors)
+        loadings : np.ndarray, shape (``n_features``, ``n_factors``)
             The loadings matrix
-        rotation_mtx : numpy array, shape (n_factors, n_factors)
+        rotation_mtx : np.ndarray, shape (``n_factors``, ``n_factors``)
             The rotation matrix
-        psi : numpy array or None, shape (n_factors, n_factors)
-            The factor correlations
-            matrix. This only exists
-            if the rotation is oblique.
+        psi : np.ndarray, shape (``n_factors``, ``n_factors``)
+            The factor correlations matrix. This only exists if the
+            rotation is oblique.
         """
-
-        if method == 'oblimin':
+        if method == "oblimin":
             objective = self._oblimin_obj
-        elif method == 'quartimin':
+        elif method == "quartimin":
             objective = self._quartimin_obj
-        elif method == 'geomin_obl':
+        elif method == "geomin_obl":
             objective = self._geomin_obj
 
         # initialize the rotation matrix
-        n_rows, n_cols = loadings.shape
+        _, n_cols = loadings.shape
         rotation_matrix = np.eye(n_cols)
 
         # default alpha level
@@ -335,25 +331,26 @@ def _oblique(self, loadings, method):
         new_loadings = np.dot(loadings, rotation_matrix_inv.T)
 
         obj = objective(new_loadings)
-        gradient = -np.dot(new_loadings.T, np.dot(obj['grad'], rotation_matrix_inv)).T
-        criterion = obj['criterion']
+        gradient = -np.dot(new_loadings.T, np.dot(obj["grad"], rotation_matrix_inv)).T
+        criterion = obj["criterion"]
 
         obj_t = objective(new_loadings)
 
         # main iteration loop, up to `max_iter`, calculate the gradient
-        for i in range(0, self.max_iter + 1):
-            gradient_new = gradient - np.dot(rotation_matrix,
-                                             np.diag(np.dot(np.ones(gradient.shape[0]),
-                                                            rotation_matrix * gradient)))
+        for _ in range(0, self.max_iter + 1):
+            gradient_new = gradient - np.dot(
+                rotation_matrix,
+                np.diag(np.dot(np.ones(gradient.shape[0]), rotation_matrix * gradient)),
+            )
             s = np.sqrt(np.sum(np.diag(np.dot(gradient_new.T, gradient_new))))
 
-            if (s < self.tol):
+            if s < self.tol:
                 break
 
             alpha = 2 * alpha
 
             # calculate the Hessian of the objective function
-            for j in range(0, 11):
+            for _ in range(0, 11):
                 X = rotation_matrix - alpha * gradient_new
 
                 v = 1 / np.sqrt(np.dot(np.ones(X.shape[0]), X**2))
@@ -361,17 +358,19 @@ def _oblique(self, loadings, method):
                 new_loadings = np.dot(loadings, np.linalg.inv(new_rotation_matrix).T)
 
                 obj_t = objective(new_loadings)
-                improvement = criterion - obj_t['criterion']
+                improvement = criterion - obj_t["criterion"]
 
-                if (improvement > 0.5 * s**2 * alpha):
+                if improvement > 0.5 * s**2 * alpha:
                     break
 
                 alpha = alpha / 2
 
             rotation_matrix = new_rotation_matrix
-            criterion = obj_t['criterion']
-            gradient = -np.dot(np.dot(new_loadings.T, obj_t['grad']),
-                               np.linalg.inv(new_rotation_matrix)).T
+            criterion = obj_t["criterion"]
+            gradient = -np.dot(
+                np.dot(new_loadings.T, obj_t["grad"]),
+                np.linalg.inv(new_rotation_matrix),
+            ).T
 
         # calculate phi
         phi = np.dot(rotation_matrix.T, rotation_matrix)
@@ -382,39 +381,38 @@ def _oblique(self, loadings, method):
 
     def _orthogonal(self, loadings, method):
         """
-        A generic function for performing
-        all orthogonal rotations, except for
-        varimax, which is implemented
-        separately.
+        Perform orthogonal rotations, except 'varimax'.
+
+        A generic function for performing all orthogonal rotations, except for
+        varimax, which is implemented separately.
 
         Parameters
         ----------
-        loadings : numpy array
+        loadings : np.ndarray
             The loading matrix
         method : str
             The orthogonal rotation method to use.
 
         Returns
         -------
-        loadings : numpy array, shape (n_features, n_factors)
+        loadings : np.ndarray
             The loadings matrix
-        rotation_mtx : numpy array, shape (n_factors, n_factors)
+        rotation_mtx : np.ndarray, shape (``n_factors``, ``n_factors``)
             The rotation matrix
         """
-
-        if method == 'oblimax':
+        if method == "oblimax":
             objective = self._oblimax_obj
-        elif method == 'quartimax':
+        elif method == "quartimax":
             objective = self._quartimax_obj
-        elif method == 'equamax':
+        elif method == "equamax":
             objective = self._equamax_obj
-        elif method == 'geomin_ort':
+        elif method == "geomin_ort":
             objective = self._geomin_obj
 
         arr = loadings.copy()
 
         # initialize the rotation matrix
-        n_rows, n_cols = arr.shape
+        _, n_cols = arr.shape
         rotation_matrix = np.eye(n_cols)
 
         # default alpha level
@@ -422,40 +420,40 @@ def _orthogonal(self, loadings, method):
         new_loadings = np.dot(arr, rotation_matrix)
 
         obj = objective(new_loadings)
-        gradient = np.dot(arr.T, obj['grad'])
-        criterion = obj['criterion']
+        gradient = np.dot(arr.T, obj["grad"])
+        criterion = obj["criterion"]
 
         obj_t = objective(new_loadings)
 
         # main iteration loop, up to `max_iter`, calculate the gradient
-        for i in range(0, self.max_iter + 1):
+        for _ in range(0, self.max_iter + 1):
             M = np.dot(rotation_matrix.T, gradient)
             S = (M + M.T) / 2
             gradient_new = gradient - np.dot(rotation_matrix, S)
             s = np.sqrt(np.sum(np.diag(np.dot(gradient_new.T, gradient_new))))
 
-            if (s < self.tol):
+            if s < self.tol:
                 break
 
             alpha = 2 * alpha
 
             # calculate the Hessian of the objective function
-            for j in range(0, 11):
+            for _ in range(0, 11):
                 X = rotation_matrix - alpha * gradient_new
-                U, D, V = np.linalg.svd(X)
+                U, _, V = np.linalg.svd(X)
                 new_rotation_matrix = np.dot(U, V)
                 new_loadings = np.dot(arr, new_rotation_matrix)
 
                 obj_t = objective(new_loadings)
 
-                if (obj_t['criterion'] < (criterion - 0.5 * s**2 * alpha)):
+                if obj_t["criterion"] < (criterion - 0.5 * s**2 * alpha):
                     break
 
                 alpha = alpha / 2
 
             rotation_matrix = new_rotation_matrix
-            criterion = obj_t['criterion']
-            gradient = np.dot(arr.T, obj_t['grad'])
+            criterion = obj_t["criterion"]
+            gradient = np.dot(arr.T, obj_t["grad"])
 
         # convert loadings matrix to data frame
         loadings = new_loadings.copy()
@@ -463,20 +461,19 @@ def _orthogonal(self, loadings, method):
 
     def _varimax(self, loadings):
         """
-        Perform varimax (orthogonal) rotation, with optional
-        Kaiser normalization.
+        Perform varimax (orthogonal) rotation, with optional Kaiser normalization.
 
         Parameters
         ----------
         loadings : array-like
-            The loading matrix
+            The loading matrix.
 
         Returns
         -------
-        loadings : numpy array, shape (n_features, n_factors)
-            The loadings matrix
-        rotation_mtx : numpy array, shape (n_factors, n_factors)
-            The rotation matrix
+        loadings : np.ndarray, shape (``n_features``, ``n_factors``)
+            The loadings matrix.
+        rotation_mtx : np.ndarray, shape (``n_factors``, ``n_factors``)
+            The rotation matrix.
         """
         X = loadings.copy()
         n_rows, n_cols = X.shape
@@ -486,7 +483,9 @@ def _varimax(self, loadings):
         # normalize the loadings matrix
         # using sqrt of the sum of squares (Kaiser)
         if self.normalize:
-            normalized_mtx = np.apply_along_axis(lambda x: np.sqrt(np.sum(x**2)), 1, X.copy())
+            normalized_mtx = np.apply_along_axis(
+                lambda x: np.sqrt(np.sum(x**2)), 1, X.copy()
+            )
             X = (X.T / normalized_mtx).T
 
         # initialize the rotation matrix
@@ -535,8 +534,7 @@ def _varimax(self, loadings):
 
     def _promax(self, loadings):
         """
-        Perform promax (oblique) rotation, with optional
-        Kaiser normalization.
+        Perform promax (oblique) rotation, with optional Kaiser normalization.
 
         Parameters
         ----------
@@ -545,17 +543,16 @@ def _promax(self, loadings):
 
         Returns
         -------
-        loadings : numpy array, shape (n_features, n_factors)
+        loadings : np.ndarray, shape (``n_features``, ``n_factors``)
             The loadings matrix
-        rotation_mtx : numpy array, shape (n_factors, n_factors)
+        rotation_mtx : np.ndarray, shape (``n_factors``, ``n_factors``)
             The rotation matrix
-        psi : numpy array or None, shape (n_factors, n_factors)
-            The factor correlations
-            matrix. This only exists
-            if the rotation is oblique.
+        psi : np.ndarray or None, shape (``n_factors``, ``n_factors``)
+            The factor correlations matrix. This only exists if the rotation
+            is oblique.
         """
         X = loadings.copy()
-        n_rows, n_cols = X.shape
+        _, n_cols = X.shape
         if n_cols < 2:
             return X
 
@@ -572,7 +569,7 @@ def _promax(self, loadings):
 
         # first get varimax rotation
         X, rotation_mtx = self._varimax(weights)
-        Y = X * np.abs(X)**(self.power - 1)
+        Y = X * np.abs(X) ** (self.power - 1)
 
         # fit linear regression model
         coef = np.dot(np.linalg.inv(np.dot(X.T, X)), np.dot(X.T, Y))
@@ -603,13 +600,13 @@ def _promax(self, loadings):
 
     def fit(self, X, y=None):
         """
-        Computes the factor rotation.
+        Compute the factor rotation.
 
         Parameters
         ----------
         X : array-like
-            The factor loading matrix (n_features, n_factors)
-        y : Ignored
+            The factor loading matrix, shape (``n_features``, ``n_factors``)
+        y : ignored
 
         Returns
         -------
@@ -630,26 +627,23 @@ def fit(self, X, y=None):
 
     def fit_transform(self, X, y=None):
         """
-        Computes the factor rotation,
-        and returns the new loading matrix.
+        Compute the factor rotation, and return the new loading matrix.
 
         Parameters
         ----------
         X : array-like
-            The factor loading matrix (n_features, n_factors)
+            The factor loading matrix, shape (``n_features``, ``n_factors``)
         y : Ignored
 
         Returns
         -------
-        loadings_ : numpy array, shape (n_features, n_factors)
-            The loadings matrix
-            (n_features, n_factors)
+        loadings_ : np,ndarray, shape (``n_features``, ``n_factors``)
+            The loadings matrix.
 
         Raises
         ------
         ValueError
-            If the `method` is not in the list of
-            acceptable methods.
+            If ``method`` is not in the list of acceptable methods.
 
         Example
         -------
@@ -676,29 +670,27 @@ def fit_transform(self, X, y=None):
         # for oblique rotations
         phi = None
         method = self.method.lower()
-        if method == 'varimax':
-            (new_loadings,
-             new_rotation_mtx) = self._varimax(X)
+        if method == "varimax":
+            (new_loadings, new_rotation_mtx) = self._varimax(X)
 
-        elif method == 'promax':
-            (new_loadings,
-             new_rotation_mtx,
-             phi) = self._promax(X)
+        elif method == "promax":
+            (new_loadings, new_rotation_mtx, phi) = self._promax(X)
 
         elif method in OBLIQUE_ROTATIONS:
-            (new_loadings,
-             new_rotation_mtx,
-             phi) = self._oblique(X, method)
+            (new_loadings, new_rotation_mtx, phi) = self._oblique(X, method)
 
         elif method in ORTHOGONAL_ROTATIONS:
-            (new_loadings,
-             new_rotation_mtx) = self._orthogonal(X, method)
+            (new_loadings, new_rotation_mtx) = self._orthogonal(X, method)
 
         else:
-            raise ValueError("The value for `method` must be one of the "
-                             "following: {}.".format(', '.join(POSSIBLE_ROTATIONS)))
-
-        (self.loadings_,
-         self.rotation_,
-         self.phi_) = new_loadings, new_rotation_mtx, phi
+            raise ValueError(
+                "The value for `method` must be one of the "
+                "following: {}.".format(", ".join(POSSIBLE_ROTATIONS))
+            )
+
+        (self.loadings_, self.rotation_, self.phi_) = (
+            new_loadings,
+            new_rotation_mtx,
+            phi,
+        )
         return self.loadings_
diff --git a/factor_analyzer/utils.py b/factor_analyzer/utils.py
index 8479e2e..afbb36d 100644
--- a/factor_analyzer/utils.py
+++ b/factor_analyzer/utils.py
@@ -1,10 +1,10 @@
 """
 Utility functions, used primarily by the confirmatory factor analysis module.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
 :author: Nitin Madnani (nmadnani@ets.org)
 :organization: Educational Testing Service
-:date: 2021-10-18
+:date: 2022-09-05
 """
 import warnings
 
@@ -14,27 +14,25 @@
 
 def inv_chol(x, logdet=False):
     """
-    Calculate inverse using cholesky.
-    Optionally, calculate the log determinant
-    of the cholesky.
+    Calculate matrix inverse using Cholesky decomposition.
+
+    Optionally, calculate the log determinant of the Cholesky.
 
     Parameters
     ----------
     x : array-like
         The matrix to invert.
     logdet : bool, optional
-        Whether to calculate the
-        log determinant, instead of
-        the inverse.
-        Defaults to False.
+        Whether to calculate the log determinant, instead of the inverse.
+        Defaults to ``False``.
 
     Returns
     -------
     chol_inv : array-like
-        The inverted matrix
+        The inverted matrix.
     chol_logdet : array-like or None
-        The log determinant, if `logdet=True`;
-        otherwise, None.
+        The log determinant, if ``logdet`` was ``True``, otherwise, ``None``.
+
     """
     chol = cholesky(x, lower=True)
 
@@ -94,57 +92,52 @@ def corr(x):
     return r
 
 
-def apply_impute_nan(x, how='mean'):
+def apply_impute_nan(x, how="mean"):
     """
-    Apply a function to impute NaN values
-    with the mean or median.
+    Apply a function to impute ``np.nan`` values with the mean or the median.
 
     Parameters
     ----------
     x : array-like
-        A 1-D array to impute.
+        The 1-D array to impute.
     how : str, optional
-        Whether to impute the 'mean'
-        or 'median'.
+        Whether to impute the 'mean' or 'median'.
         Defaults to 'mean'.
 
     Returns
     -------
-    x : numpy array
-        The array, with missing values imputed.
+    x : np.ndarray
+        The array, with the missing values imputed.
     """
-
-    if how == 'mean':
+    if how == "mean":
         x[np.isnan(x)] = np.nanmean(x)
-    elif how == 'median':
+    elif how == "median":
         x[np.isnan(x)] = np.nanmedian(x)
     return x
 
 
-def impute_values(x, how='mean'):
+def impute_values(x, how="mean"):
     """
-    Impute NaN values with the mean or median,
-    or drop rows with NaN values.
+    Impute ``np.nan`` values with the mean or median, or drop the containing rows.
 
     Parameters
     ----------
     x : array-like
         An array to impute.
     how : str, optional
-        Whether to impute the 'mean'
-        or 'median'.
+        Whether to impute the 'mean' or 'median'.
         Defaults to 'mean'.
 
     Returns
     -------
-    x : numpy array
-        The array, with missing values imputed or dropped.
+    x : np.ndarray
+        The array, with the missing values imputed or with rows dropped.
     """
     # impute mean or median, if `how` is set to 'mean' or 'median'
-    if how in ['mean', 'median']:
+    if how in ["mean", "median"]:
         x = np.apply_along_axis(apply_impute_nan, 0, x, how=how)
     # drop missing if `how` is set to 'drop'
-    elif how == 'drop':
+    elif how == "drop":
         x = x[~np.isnan(x).any(1), :].copy()
     return x
 
@@ -152,24 +145,23 @@ def impute_values(x, how='mean'):
 def smc(corr_mtx, sort=False):
     """
     Calculate the squared multiple correlations.
-    This is equivalent to regressing each variable
-    on all others and calculating the r-squared values.
+
+    This is equivalent to regressing each variable on all others and
+    calculating the r-squared values.
 
     Parameters
     ----------
     corr_mtx : array-like
         The correlation matrix used to calculate SMC.
     sort : bool, optional
-        Whether to sort the values for SMC
-        before returning.
-        Defaults to False.
+        Whether to sort the values for SMC before returning.
+        Defaults to ``False``.
 
     Returns
     -------
-    smc : numpy array
+    smc : np.ndarray
         The squared multiple correlations matrix.
     """
-
     corr_inv = np.linalg.inv(corr_mtx)
     smc = 1 - 1 / np.diag(corr_inv)
 
@@ -180,7 +172,9 @@ def smc(corr_mtx, sort=False):
 
 def covariance_to_correlation(m):
     """
-    This is a port of the R `cov2cor` function.
+    Compute cross-correlations from the given covariance matrix.
+
+    This is a port of R ``cov2cor()`` function.
 
     Parameters
     ----------
@@ -189,7 +183,7 @@ def covariance_to_correlation(m):
 
     Returns
     -------
-    retval : numpy array
+    retval : np.ndarray
         The cross-correlation matrix.
 
     Raises
@@ -197,11 +191,10 @@ def covariance_to_correlation(m):
     ValueError
         If the input matrix is not square.
     """
-
     # make sure the matrix is square
     numrows, numcols = m.shape
     if not numrows == numcols:
-        raise ValueError('Input matrix must be square')
+        raise ValueError("Input matrix must be square")
 
     Is = np.sqrt(1 / np.diag(m))
     retval = Is * m * np.repeat(Is, numrows).reshape(numrows, numrows)
@@ -211,9 +204,11 @@ def covariance_to_correlation(m):
 
 def partial_correlations(x):
     """
-    This is a python port of the `pcor` function implemented in
-    the `ppcor` R package, which computes partial correlations
-    of each pair of variables in the given array, excluding all
+    Compute partial correlations between variable pairs.
+
+    This is a python port of the ``pcor()`` function implemented in
+    the ``ppcor`` R package, which computes partial correlations
+    for each pair of variables in the given array, excluding all
     other variables.
 
     Parameters
@@ -223,7 +218,7 @@ def partial_correlations(x):
 
     Returns
     -------
-    pcor : numpy array
+    pcor : np.ndarray
         An array containing the partial correlations of of each
         pair of variables in the given array, excluding all other
         variables.
@@ -246,10 +241,12 @@ def partial_correlations(x):
             icvx = np.linalg.inv(x_cov)
         except AssertionError:
             icvx = np.linalg.pinv(x_cov)
-            warnings.warn('The inverse of the variance-covariance matrix '
-                          'was calculated using the Moore-Penrose generalized '
-                          'matrix inversion, due to its determinant being at '
-                          'or very close to zero.')
+            warnings.warn(
+                "The inverse of the variance-covariance matrix "
+                "was calculated using the Moore-Penrose generalized "
+                "matrix inversion, due to its determinant being at "
+                "or very close to zero."
+            )
         except np.linalg.LinAlgError:
             icvx = empty_array
 
@@ -260,9 +257,7 @@ def partial_correlations(x):
 
 def unique_elements(seq):
     """
-    Get the first unique instance of every
-    element of the list, while maintaining
-    the original order of those instances.
+    Get first unique instance of every list element, while maintaining order.
 
     Parameters
     ----------
@@ -280,23 +275,21 @@ def unique_elements(seq):
 
 def fill_lower_diag(x):
     """
-    Fill the lower diagonal of a square matrix,
-    given a 1d input array.
+    Fill the lower diagonal of a square matrix, given a 1-D input array.
 
     Parameters
     ----------
     x : array-like
-        The flattened input matrix that will be used to fill
-        the lower diagonal of the square matrix.
+        The flattened input matrix that will be used to fill the lower
+        diagonal of the square matrix.
 
     Returns
     -------
-    out : numpy array
-        The output square matrix, with the lower
-        diagonal filled by x.
+    out : np.ndarray
+        The output square matrix, with the lower diagonal filled by x.
 
-    Reference
-    ---------
+    References
+    ----------
     [1] https://stackoverflow.com/questions/51439271/
         convert-1d-array-to-lower-triangular-matrix
     """
@@ -310,30 +303,29 @@ def fill_lower_diag(x):
 
 def merge_variance_covariance(variances, covariances=None):
     """
-    Merge the variance and covariances into a single
-    variance-covariance matrix.
+    Merge variances and covariances into a single variance-covariance matrix.
 
     Parameters
     ----------
     variances : array-like
-        The variances that will be used to fill the diagonal
-        of the square matrix.
+        The variances that will be used to fill the diagonal of the
+        square matrix.
     covariances : array-like or None, optional
-        The flattened input matrix that will be used to fill
-        the lower and upper diagonal of the square matrix. If
-        None, then only the variances will be used.
-        Defaults to None.
+        The flattened input matrix that will be used to fill the lower and
+        upper diagonal of the square matrix. If None, then only the variances
+        will be used.
+        Defaults to ``None``.
 
     Returns
     -------
-    variance_covariance : numpy array
+    variance_covariance : np.ndarray
         The variance-covariance matrix.
     """
-    variances = (variances if len(variances.shape) == 1
-                 else np.squeeze(variances, axis=1))
+    variances = (
+        variances if len(variances.shape) == 1 else np.squeeze(variances, axis=1)
+    )
     if covariances is None:
-        variance_covariance = np.zeros((variances.shape[0],
-                                        variances.shape[0]))
+        variance_covariance = np.zeros((variances.shape[0], variances.shape[0]))
     else:
         variance_covariance = fill_lower_diag(covariances)
         variance_covariance += variance_covariance.T
@@ -343,7 +335,7 @@ def merge_variance_covariance(variances, covariances=None):
 
 def get_first_idxs_from_values(x, eq=1, use_columns=True):
     """
-    Get the fixed index
+    Get the indexes  for a given value.
 
     Parameters
     ----------
@@ -353,10 +345,9 @@ def get_first_idxs_from_values(x, eq=1, use_columns=True):
         The given value to find.
         Defaults to 1.
     use_columns : bool, optional
-        Whether to get the first indexes using
-        The columns. If False, then use the rows
-        instead.
-        Defaults to True
+        Whether to get the first indexes using the columns.
+        If ``False``, then use the rows instead.
+        Defaults to ``True``.
 
     Returns
     -------
@@ -379,31 +370,31 @@ def get_first_idxs_from_values(x, eq=1, use_columns=True):
 
 def get_free_parameter_idxs(x, eq=1):
     """
-    Get the free parameter indexes from
-    the flattened matrix.
+    Get the free parameter indices from the flattened matrix.
 
     Parameters
     ----------
     x : array-like
         The input matrix.
     eq : str or int, optional
-        The value that free parameters
-        should be equal to. NaN fields
+        The value that free parameters should be equal to. ``np.nan`` fields
         will be populated with this value.
-        Defaults to 1
+        Defaults to 1.
 
     Returns
     -------
-    idx : numpy array
+    idx : np.ndarray
         The free parameter indexes.
     """
     x[np.isnan(x)] = eq
-    x = x.flatten(order='F')
+    x = x.flatten(order="F")
     return np.where(x == eq)[0]
 
 
 def duplication_matrix(n=1):
     """
+    Calculate the duplication matrix.
+
     A function to create the duplication matrix (Dn), which is
     the unique n2 × n(n+1)/2 matrix which, for any n × n symmetric
     matrix A, transforms vech(A) into vec(A), as in Dn vech(A) = vec(A).
@@ -416,27 +407,28 @@ def duplication_matrix(n=1):
 
     Returns
     -------
-    duplication_matrix : numpy array
+    duplication_matrix : np.ndarray
         The duplication matrix.
 
     Raises`
     ------
     ValueError
-        If `n` is not a positive integer greater than 1.
+        If ``n`` is not a positive integer greater than 1.
 
     References
     ----------
     https://en.wikipedia.org/wiki/Duplication_and_elimination_matrices
     """
     if n < 1:
-        raise ValueError('The argument `n` must be a '
-                         'positive integer greater than 1.')
+        raise ValueError(
+            "The argument `n` must be a " "positive integer greater than 1."
+        )
 
     dup = np.zeros((int(n * n), int(n * (n + 1) / 2)))
     count = 0
     for j in range(n):
         dup[j * n + j, count + j] = 1
-        if (j < n - 1):
+        if j < n - 1:
             for i in range(j + 1, n):
                 dup[j * n + i, count + i] = 1
                 dup[i * n + j, count + i] = 1
@@ -446,25 +438,23 @@ def duplication_matrix(n=1):
 
 def duplication_matrix_pre_post(x):
     """
-    Given an input symmetric matrix,
-    transform perform the pre-post duplication.
+    Transform given input symmetric matrix using pre-post duplication.
 
     Parameters
     ----------
     x : array-like
-        The input matrix
+        The input matrix.
 
     Returns
     -------
-    out : numpy array
+    out : np.ndarray
         The transformed matrix.
 
     Raises
     ------
     AssertionError
-        If `x` is not symmetric.
+        If ``x`` is not symmetric.
     """
-
     assert x.shape[0] == x.shape[1]
 
     n2 = x.shape[1]
@@ -483,9 +473,10 @@ def duplication_matrix_pre_post(x):
 
 def commutation_matrix(p, q):
     """
-    Calculate the commutation matrix, which transforms
-    the vectorized form of the matrix into the vectorized
-    form of its transpose.
+    Calculate the commutation matrix.
+
+    This matrix transforms the vectorized form of the matrix into the
+    vectorized form of its transpose.
 
     Parameters
     ----------
@@ -496,7 +487,7 @@ def commutation_matrix(p, q):
 
     Returns
     -------
-    commutation_matrix : numpy array
+    commutation_matrix : np.ndarray
         The commutation matrix
 
     References
@@ -504,14 +495,13 @@ def commutation_matrix(p, q):
     https://en.wikipedia.org/wiki/Commutation_matrix
     """
     identity = np.eye(p * q)
-    indices = np.arange(p * q).reshape((p, q), order='F')
+    indices = np.arange(p * q).reshape((p, q), order="F")
     return identity.take(indices.ravel(), axis=0)
 
 
 def get_symmetric_lower_idxs(n=1, diag=True):
     """
-    Get the indexes for the lower triangle of
-    a symmetric matrix.
+    Get the indices for the lower triangle of a symmetric matrix.
 
     Parameters
     ----------
@@ -523,8 +513,8 @@ def get_symmetric_lower_idxs(n=1, diag=True):
 
     Returns
     -------
-    indexes : numpy array
-        The indexes for the lower triangle.
+    indices : np.ndarray
+        The indices for the lower triangle.
     """
     rows = np.repeat(np.arange(n), n).reshape(n, n)
     cols = rows.T
@@ -535,8 +525,7 @@ def get_symmetric_lower_idxs(n=1, diag=True):
 
 def get_symmetric_upper_idxs(n=1, diag=True):
     """
-    Get the indexes for the upper triangle of
-    a symmetric matrix.
+    Get the indices for the upper triangle of a symmetric matrix.
 
     Parameters
     ----------
@@ -548,8 +537,8 @@ def get_symmetric_upper_idxs(n=1, diag=True):
 
     Returns
     -------
-    indexes : numpy array
-        The indexes for the upper triangle.
+    indices : np.ndarray
+        The indices for the upper triangle.
     """
     rows = np.repeat(np.arange(n), n).reshape(n, n)
     cols = rows.T
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..b2004ea
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,5 @@
+[tool.isort]
+profile = "black"
+
+[pydocstyle]
+convention = "pep257"
diff --git a/requirements.txt b/requirements.txt
index 1a9357d..38a6827 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,5 @@
 pandas
 scipy
 numpy
-scikit-learn
\ No newline at end of file
+scikit-learn
+pre-commit
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..2f394d9
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,7 @@
+[wheel]
+universal = 1
+
+[flake8]
+max_line_length = 101
+exclude =
+    factor_analyzer/__init__.py
diff --git a/setup.py b/setup.py
index 60fcf3c..363e598 100644
--- a/setup.py
+++ b/setup.py
@@ -1,42 +1,49 @@
 #!/usr/bin/env python
+"""Set up the package."""
 
 from setuptools import find_packages, setup
 
 
-def readme():
-    with open('README.rst') as f:
+def readme():  # noqa: D103
+    with open("README.rst") as f:
         return f.read()
 
 
-def requirements():
-    req_path = 'requirements.txt'
+def requirements():  # noqa: D103
+    req_path = "requirements.txt"
     with open(req_path) as f:
         reqs = f.read().splitlines()
     return reqs
 
 
-setup(name='factor_analyzer',
-      version='0.4.0',
-      description='A Factor Analysis class',
-      long_description=readme(),
-      keywords='factor analysis',
-      packages=find_packages(),
-      author="Jeremy Biggs",
-      author_email="jbiggs@ets.org",
-      url="https://github.com/EducationalTestingService/factor_analyzer",
-      install_requires=requirements(),
-      include_package_data=True,
-      classifiers=['Intended Audience :: Science/Research',
-                   'Intended Audience :: Developers',
-                   'Programming Language :: Python',
-                   'Topic :: Scientific/Engineering',
-                   'Operating System :: Microsoft :: Windows',
-                   'Operating System :: POSIX',
-                   'Operating System :: Unix',
-                   'Operating System :: MacOS',
-                   'Programming Language :: Python :: 3',
-                   'Programming Language :: Python :: 3.7',
-                   'Programming Language :: Python :: 3.8',
-                   'Programming Language :: Python :: 3.9',
-                   ],
-      zip_safe=True)
+setup(
+    name="factor_analyzer",
+    version="0.4.0",
+    description="A Factor Analysis tool written in Python",
+    long_description=readme(),
+    keywords="factor analysis",
+    packages=find_packages(),
+    author="Jeremy Biggs",
+    author_email="jeremy.m.biggs@gmail.com",
+    maintainer="Nitin Madnani",
+    maintainer_email="nmadnani@ets.org",
+    url="https://github.com/EducationalTestingService/factor_analyzer",
+    install_requires=requirements(),
+    include_package_data=True,
+    classifiers=[
+        "Intended Audience :: Science/Research",
+        "Intended Audience :: Developers",
+        "Programming Language :: Python",
+        "Topic :: Scientific/Engineering",
+        "Operating System :: Microsoft :: Windows",
+        "Operating System :: POSIX",
+        "Operating System :: Unix",
+        "Operating System :: MacOS",
+        "Programming Language :: Python :: 3",
+        "Programming Language :: Python :: 3.7",
+        "Programming Language :: Python :: 3.8",
+        "Programming Language :: Python :: 3.9",
+        "Programming Language :: Python :: 3.10",
+    ],
+    zip_safe=True,
+)


Test patch:
diff --git a/factor_analyzer/test_utils.py b/factor_analyzer/test_utils.py
index bc5ce13..0254afb 100644
--- a/factor_analyzer/test_utils.py
+++ b/factor_analyzer/test_utils.py
@@ -1,10 +1,10 @@
 """
 Utility functions used for testing.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
 :author: Nitin Madnani (nmadnani@ets.org)
 :organization: Educational Testing Service
-:date: 2021-10-18
+:date: 2022-09-05
 """
 
 import json
@@ -14,63 +14,73 @@
 
 import numpy as np
 import pandas as pd
-from factor_analyzer import ConfirmatoryFactorAnalyzer, FactorAnalyzer, ModelSpecificationParser, Rotator
+
+from factor_analyzer import (
+    ConfirmatoryFactorAnalyzer,
+    FactorAnalyzer,
+    ModelSpecificationParser,
+    Rotator,
+)
 from factor_analyzer.utils import unique_elements
 
-DATA_DIR = os.path.join('tests', 'data')
-JSON_DIR = os.path.join('tests', 'model')
-EXPECTED_DIR = os.path.join('tests', 'expected')
-
-OUTPUT_TYPES = ['value',
-                'evalues',
-                'loading',
-                'uniquenesses',
-                'communalities',
-                'structure',
-                'scores']
-
-
-def calculate_py_output(test_name,
-                        factors,
-                        method,
-                        rotation,
-                        svd_method='randomized',
-                        use_corr_matrix=False,
-                        top_dir=None):
+DATA_DIR = os.path.join("tests", "data")
+JSON_DIR = os.path.join("tests", "model")
+EXPECTED_DIR = os.path.join("tests", "expected")
+
+OUTPUT_TYPES = [
+    "value",
+    "evalues",
+    "loading",
+    "uniquenesses",
+    "communalities",
+    "structure",
+    "scores",
+]
+
+
+def calculate_py_output(
+    test_name,
+    factors,
+    method,
+    rotation,
+    svd_method="randomized",
+    use_corr_matrix=False,
+    top_dir=None,
+):
     """
-    Use the `FactorAnalyzer()` class to perform the factor analysis
-    and return a dictionary with relevant output for given scenario.
+    Use :class:`FactorAnalyzer()` to perform factor analysis and return results.
+
+    The results are returned in a dictionary.
 
     Parameters
     ----------
     test_name : str
-        The name of the test
+        The name of the test.
     factors : int
-        The number of factors
+        The number of factors.
     method : str
-        The rotation method
+        The rotation method.
     rotation : str
-        The type of rotation
+        The type of rotation.
     svd_method : str, optional
-        The SVD method to use
+        The SVD method to use.
         Defaults to 'randomized'
     use_corr_matrix : bool, optional
         Whether to use the correlation matrix.
-        Defaults to False.
+        Defaults to ``False``.
     top_dir : str, optional
-        The top directory for test data
-        Defaults to `DATA_DIR``
+        The top directory for test data.
+        Defaults to ``DATA_DIR``
 
     Returns
     -------
     output : dict
-        A dictionary containing the outputs
-        for all `OUTPUT_TYPES`.
+        A dictionary containing the outputs for all ``OUTPUT_TYPES``.
     """
     if top_dir is None:
         top_dir = DATA_DIR
 
-    filename = join(top_dir, test_name + '.csv')
+    filename = join(top_dir, test_name + ".csv")
     data = pd.read_csv(filename)
 
     if use_corr_matrix:
@@ -78,31 +88,34 @@ def calculate_py_output(test_name,
     else:
         X = data.copy()
 
-    rotation = None if rotation == 'none' else rotation
-    method = {'uls': 'minres'}.get(method, method)
+    rotation = None if rotation == "none" else rotation
+    method = {"uls": "minres"}.get(method, method)
 
-    fa = FactorAnalyzer(n_factors=factors, method=method, svd_method=svd_method,
-                        rotation=rotation, is_corr_matrix=use_corr_matrix)
+    fa = FactorAnalyzer(
+        n_factors=factors,
+        method=method,
+        svd_method=svd_method,
+        rotation=rotation,
+        is_corr_matrix=use_corr_matrix,
+    )
     fa.fit(X)
 
     evalues, values = fa.get_eigenvalues()
 
-    return {'value': values,
-            'evalues': evalues,
-            'structure': fa.structure_,
-            'loading': fa.loadings_,
-            'uniquenesses': fa.get_uniquenesses(),
-            'communalities': fa.get_communalities(),
-            'scores': fa.transform(data)}
-
-
-def collect_r_output(test_name,
-                     factors,
-                     method,
-                     rotation,
-                     output_types=None,
-                     top_dir=None,
-                     **kwargs):
+    return {
+        "value": values,
+        "evalues": evalues,
+        "structure": fa.structure_,
+        "loading": fa.loadings_,
+        "uniquenesses": fa.get_uniquenesses(),
+        "communalities": fa.get_communalities(),
+        "scores": fa.transform(data),
+    }
+
+
+def collect_r_output(
+    test_name, factors, method, rotation, output_types=None, top_dir=None, **kwargs
+):
     """
     Get the R output for the given scenario.
 
@@ -142,11 +155,7 @@ def collect_r_output(test_name,
     output = {}
     for output_type in output_types:
 
-        filename = '{}_{}_{}_{}_{}.csv'.format(output_type,
-                                               method,
-                                               rotation,
-                                               factors,
-                                               test_name)
+        filename = f"{output_type}_{method}_{rotation}_{factors}_{test_name}.csv"
 
         filename = os.path.join(top_dir, test_name, filename)
 
@@ -158,31 +167,29 @@ def collect_r_output(test_name,
 
 def normalize(data, absolute=False):
     """
-    Normalize the data to ensure that Python
-    and R output match. This involves ensuring
-    the headers are named consistently and
-    columns are sorted properly.
+    Normalize the data to ensure that Python and R output match.
+
+    This involves ensuring the headers are named consistently and columns
+    are sorted properly.
 
     Parameters
     ----------
     data : pd.DataFrame
         The data frame to normalize.
     absolute : bool, optional
-        Whether to take the absolute value of
-        all elements in the data frame.
-        Defaults to True
+        Whether to take the absolute value of all elements in the data frame.
+        Defaults to ``True``.
 
     Returns
     -------
     data : pd.DataFrame
         The normalized data frame.
     """
-
     # # get numeric columns, in case we are taking absolute value
     data = data.copy()
 
     if isinstance(data, pd.DataFrame):
-        data = data[[col for col in data if 'unnamed' not in col.lower()]].copy()
+        data = data[[col for col in data if "unnamed" not in col.lower()]].copy()
         data = data.select_dtypes(include=[np.number])
         data = data.values
     else:
@@ -196,11 +203,13 @@ def normalize(data, absolute=False):
     return data
 
 
-def check_close(data1, data2, rel_tol=0.0, abs_tol=0.1,
-                with_normalize=True, absolute=False):
+def check_close(
+    data1, data2, rel_tol=0.0, abs_tol=0.1, with_normalize=True, absolute=False
+):
     """
-    Check to make sure all values in two data frames
-    are close. Returns the proportion that match.
+    Check that all values in the two data frames are close.
+
+    Returns the proportion that match.
 
     Parameters
     ----------
@@ -215,9 +224,8 @@ def check_close(data1, data2, rel_tol=0.0, abs_tol=0.1,
         The absolute tolerance.
         Defaults to 0.1.
     absolute : bool, optional
-        Whether to take the absolute value of
-        all elements in the data frame.
-        Defaults to False
+        Whether to take the absolute value of all elements in the data frame.
+        Defaults to ``False``.
 
     Returns
     -------
@@ -231,38 +239,39 @@ def check_close(data1, data2, rel_tol=0.0, abs_tol=0.1,
     print(data1)
     print()
     print(data2)
-    print('------')
+    print("------")
 
-    err_msg = 'r - py: {} != {}'
+    err_msg = "r - py: {} != {}"
     assert data1.shape == data2.shape, err_msg.format(data1.shape, data2.shape)
 
     arr = np.empty(shape=data1.shape, dtype=bool)
     for i in range(data1.shape[0]):
         for j in range(data2.shape[1]):
-            check = math.isclose(data1[i, j],
-                                 data2[i, j],
-                                 rel_tol=rel_tol,
-                                 abs_tol=abs_tol)
+            check = math.isclose(
+                data1[i, j], data2[i, j], rel_tol=rel_tol, abs_tol=abs_tol
+            )
             arr[i, j] = check
 
     check = np.sum(arr) / arr.size
     return check
 
 
-def check_scenario(test_name,
-                   factors,
-                   method,
-                   rotation,
-                   ignore_value=False,
-                   ignore_communalities=False,
-                   check_scores=False,
-                   check_structure=False,
-                   use_corr_matrix=False,
-                   svd_method='randomized',
-                   data_dir=None,
-                   expected_dir=None,
-                   rel_tol=0,
-                   abs_tol=0.1):
+def check_scenario(
+    test_name,
+    factors,
+    method,
+    rotation,
+    ignore_value=False,
+    ignore_communalities=False,
+    check_scores=False,
+    check_structure=False,
+    use_corr_matrix=False,
+    svd_method="randomized",
+    data_dir=None,
+    expected_dir=None,
+    rel_tol=0,
+    abs_tol=0.1,
+):
     """
     Check all results for given scenario.
 
@@ -277,28 +286,27 @@ def check_scenario(test_name,
     rotation : str
         The type of rotation (e.g. 'varimax')
     ignore_values : bool, optional
-        Whether to ignore the `value` output type.
-        Defaults to False.
+        Whether to ignore the "value" output type.
+        Defaults to ``False``.
     ignore_communalities : bool, optional
-        Whether to ignore the `communalities` output type.
-        Defaults to False.
+        Whether to ignore the "communalities" output type.
+        Defaults to ``False``.
     check_scores : bool, optional
-        Check the factor scores
-        Defaults to False.
+        Check the factor scores.
+        Defaults to ``False``.
     check_structure : bool, optional
-        Check the structure matrix.
-        This should only be used with
-        oblique rotations.
-        Defaults to False.
+        Check the structure matrix. This should only be used with oblique
+        rotations.
+        Defaults to ``False``.
     use_corr_matrix : bool, optional
         Whether to use the correlation matrix.
-        Defaults to False.
+        Defaults to ``False``.
     data_dir : str, optional
         The directory with input data files.
-        Defaults to `DATA_DIR`.
+        Defaults to ``DATA_DIR``.
     expected_dir : str, optional
         The directory with output files.
-        Defaults to `EXPECTED_DIR`.
+        Defaults to ``EXPECTED_DIR``.
     rel_tol : float, optional
         The relative tolerance.
         Defaults to 0.0.
@@ -309,28 +317,28 @@ def check_scenario(test_name,
     Yields
     ------
     check : float
-        The proportion that match between
-        the calculated and expected.
+        The proportion that match between the calculated and expected.
     """
-
-    output_types = ['loading', 'evalues']
+    output_types = ["loading", "evalues"]
 
     if not ignore_value:
-        output_types.append('value')
+        output_types.append("value")
 
     if not ignore_communalities:
-        output_types.extend(['uniquenesses', 'communalities'])
+        output_types.extend(["uniquenesses", "communalities"])
 
     if check_scores:
-        output_types.append('scores')
+        output_types.append("scores")
 
     if check_structure:
-        output_types.append('structure')
+        output_types.append("structure")
 
-    r_output = collect_r_output(test_name, factors, method, rotation,
-                                output_types, expected_dir)
-    py_output = calculate_py_output(test_name, factors, method, rotation, svd_method,
-                                    use_corr_matrix, data_dir)
+    r_output = collect_r_output(
+        test_name, factors, method, rotation, output_types, expected_dir
+    )
+    py_output = calculate_py_output(
+        test_name, factors, method, rotation, svd_method, use_corr_matrix, data_dir
+    )
 
     for output_type in output_types:
 
@@ -340,13 +348,9 @@ def check_scenario(test_name,
         yield check_close(data1, data2, rel_tol, abs_tol)
 
 
-def check_rotation(test_name,
-                   factors,
-                   method,
-                   rotation,
-                   rel_tol=0,
-                   abs_tol=0.1,
-                   **kwargs):
+def check_rotation(
+    test_name, factors, method, rotation, rel_tol=0, abs_tol=0.1, **kwargs
+):
     """
     Check the rotation results.
 
@@ -370,21 +374,21 @@ def check_rotation(test_name,
     Returns
     ------
     check : float
-        The proportion that match between
-        the calculated and expected.
+        The proportion that match between the calculated and expected.
     """
-
-    r_input = collect_r_output(test_name, factors, method,
-                               'none', output_types=['loading'])
-    r_loading = r_input['loading']
+    r_input = collect_r_output(
+        test_name, factors, method, "none", output_types=["loading"]
+    )
+    r_loading = r_input["loading"]
     r_loading = normalize(r_loading, absolute=False)
 
     rotator = Rotator(method=rotation, **kwargs)
     rotated_loading = rotator.fit_transform(r_loading)
 
-    r_output = collect_r_output(test_name, factors, method, rotation,
-                                output_types=['loading'])
-    expected_loading = r_output['loading']
+    r_output = collect_r_output(
+        test_name, factors, method, rotation, output_types=["loading"]
+    )
+    expected_loading = r_output["loading"]
 
     data1 = normalize(rotated_loading)
     data2 = normalize(expected_loading)
@@ -392,21 +396,23 @@ def check_rotation(test_name,
     return check_close(data1, data2, rel_tol, abs_tol)
 
 
-def calculate_py_output_cfa(json_name,
-                            data_name,
-                            is_cov=False,
-                            data_dir=None,
-                            json_dir=None,
-                            n_obs=None,
-                            **kwargs):
-
+def calculate_py_output_cfa(
+    json_name,
+    data_name,
+    is_cov=False,
+    data_dir=None,
+    json_dir=None,
+    n_obs=None,
+    **kwargs,
+):
+    """Use :class:`ConfirmatoryFactorAnalyzer()` to perform factor analysis and return results."""
     if data_dir is None:
         data_dir = DATA_DIR
     if json_dir is None:
         json_dir = JSON_DIR
 
-    filename = join(data_dir, data_name + '.csv')
-    jsonname = join(json_dir, json_name + '.json')
+    filename = join(data_dir, data_name + ".csv")
+    jsonname = join(json_dir, json_name + ".json")
     data = pd.read_csv(filename, **kwargs)
 
     if n_obs is None and not is_cov:
@@ -417,62 +423,75 @@ def calculate_py_output_cfa(json_name,
 
     columns = unique_elements([v for f in model.values() for v in f])
     data = data[columns].copy()
-    model_spec = ModelSpecificationParser.parse_model_specification_from_dict(data, model)
+    model_spec = ModelSpecificationParser.parse_model_specification_from_dict(
+        data, model
+    )
 
-    cfa = ConfirmatoryFactorAnalyzer(model_spec,
-                                     n_obs=n_obs,
-                                     is_cov_matrix=is_cov,
-                                     disp=False)
+    cfa = ConfirmatoryFactorAnalyzer(
+        model_spec, n_obs=n_obs, is_cov_matrix=is_cov, disp=False
+    )
     cfa.fit(data.values)
     transform = cfa.transform(data.values)
 
-    (loadingsse,
-     errorcovsse) = cfa.get_standard_errors()
+    (loadingsse, errorcovsse) = cfa.get_standard_errors()
 
-    outputs = {'loadings': cfa.loadings_.copy(),
-               'errorvars': cfa.error_vars_.copy(),
-               'factorcovs': cfa.factor_varcovs_.copy(),
-               'loadingsse': loadingsse.copy(),
-               'errorvarsse': errorcovsse.copy(),
-               'transform': transform.copy()}
+    outputs = {
+        "loadings": cfa.loadings_.copy(),
+        "errorvars": cfa.error_vars_.copy(),
+        "factorcovs": cfa.factor_varcovs_.copy(),
+        "loadingsse": loadingsse.copy(),
+        "errorvarsse": errorcovsse.copy(),
+        "transform": transform.copy(),
+    }
 
     return outputs, cfa.model.n_factors
 
 
-def check_cfa(json_name_input,
-              data_name_input,
-              data_name_expected=None,
-              is_cov=False,
-              n_obs=None,
-              rel_tol=0,
-              abs_tol=0,
-              data_dir=None,
-              json_dir=None,
-              expected_dir=None,
-              **kwargs):
-
+def check_cfa(
+    json_name_input,
+    data_name_input,
+    data_name_expected=None,
+    is_cov=False,
+    n_obs=None,
+    rel_tol=0,
+    abs_tol=0,
+    data_dir=None,
+    json_dir=None,
+    expected_dir=None,
+    **kwargs,
+):
+    """Check the Python CFA output against R."""
     if data_name_expected is None:
         data_name_expected = json_name_input
 
-    output_types = ['loadings', 'errorvars', 'factorcovs',
-                    'loadingsse', 'errorvarsse', 'transform']
-
-    (outputs_p,
-     factors) = calculate_py_output_cfa(json_name_input,
-                                        data_name_input,
-                                        is_cov=is_cov,
-                                        n_obs=n_obs,
-                                        data_dir=data_dir,
-                                        json_dir=json_dir,
-                                        **kwargs)
-
-    outputs_r = collect_r_output(data_name_expected,
-                                 factors,
-                                 'cfa',
-                                 'none',
-                                 output_types=output_types,
-                                 top_dir=expected_dir,
-                                 header=None)
+    output_types = [
+        "loadings",
+        "errorvars",
+        "factorcovs",
+        "loadingsse",
+        "errorvarsse",
+        "transform",
+    ]
+
+    (outputs_p, factors) = calculate_py_output_cfa(
+        json_name_input,
+        data_name_input,
+        is_cov=is_cov,
+        n_obs=n_obs,
+        data_dir=data_dir,
+        json_dir=json_dir,
+        **kwargs,
+    )
+
+    outputs_r = collect_r_output(
+        data_name_expected,
+        factors,
+        "cfa",
+        "none",
+        output_types=output_types,
+        top_dir=expected_dir,
+        header=None,
+    )
 
     for output_type in output_types:
 
@@ -481,7 +500,6 @@ def check_cfa(json_name_input,
         print(data1)
         print(data2)
 
-        yield check_close(data1, data2,
-                          rel_tol=rel_tol,
-                          abs_tol=abs_tol,
-                          with_normalize=True)
+        yield check_close(
+            data1, data2, rel_tol=rel_tol, abs_tol=abs_tol, with_normalize=True
+        )
diff --git a/tests/generate_r_output.r b/tests/generate_r_output.r
index 20a97c3..17d0003 100644
--- a/tests/generate_r_output.r
+++ b/tests/generate_r_output.r
@@ -78,10 +78,10 @@ for (n in args$n_factors) {
           evalues <- res$e.values;
           uniquenesses <- res$uniquenesses;
           communalities <- res$communalities;
-          
+
           info <- list('loading' = loadings,
                        'value' = values,
-                       'evalues' = evalues, 
+                       'evalues' = evalues,
                        'uniquenesses' = uniquenesses,
                        'communalities' = communalities)
           for (name in names(info)) {
@@ -97,4 +97,4 @@ for (n in args$n_factors) {
           }
         }
     }
-}
\ No newline at end of file
+}
diff --git a/tests/model/test12.json b/tests/model/test12.json
index edd85be..39d90e4 100644
--- a/tests/model/test12.json
+++ b/tests/model/test12.json
@@ -1,2 +1,2 @@
 {"Verbal": ["english", "vocab", "socsci"],
- "Quant": ["socsci", "math", "natsci"]}
\ No newline at end of file
+ "Quant": ["socsci", "math", "natsci"]}
diff --git a/tests/test_expected_confirmatory_factor_analyzer.py b/tests/test_expected_confirmatory_factor_analyzer.py
index b13e8f0..fb9b7f4 100644
--- a/tests/test_expected_confirmatory_factor_analyzer.py
+++ b/tests/test_expected_confirmatory_factor_analyzer.py
@@ -1,82 +1,80 @@
 """
-Tests for ConfirmatoryFactorAnalyzer class
+Tests for ``ConfirmatoryFactorAnalyzer`` class.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
-:date: 02/12/2019
-:organization: ETS
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
+:author: Nitin Madnani (nmadnani@ets.org)
+:organization: Educational Testing Service
+:date: 2022-09-05
 """
 import numpy as np
 import pandas as pd
-from factor_analyzer.test_utils import check_cfa
-from factor_analyzer.confirmatory_factor_analyzer import (ModelSpecification,
-                                                          ModelSpecificationParser,
-                                                          ConfirmatoryFactorAnalyzer)
 from nose.tools import eq_, raises
 from numpy.testing import assert_array_equal
 
+from factor_analyzer.confirmatory_factor_analyzer import (
+    ConfirmatoryFactorAnalyzer,
+    ModelSpecification,
+    ModelSpecificationParser,
+)
+from factor_analyzer.test_utils import check_cfa
 
 THRESHOLD = 1.0
 
 
-def test_11_cfa():
+def test_11_cfa():  # noqa: D103
 
-    json_name_input = 'test11'
-    data_name_input = 'test11'
+    json_name_input = "test11"
+    data_name_input = "test11"
 
-    for check in check_cfa(json_name_input,
-                           data_name_input,
-                           rel_tol=0.1):
+    for check in check_cfa(json_name_input, data_name_input, rel_tol=0.1):
         assert check >= THRESHOLD
 
 
-def test_12_cfa():
+def test_12_cfa():  # noqa: D103
 
-    json_name_input = 'test12'
-    data_name_input = 'test12'
+    json_name_input = "test12"
+    data_name_input = "test12"
 
-    for check in check_cfa(json_name_input,
-                           data_name_input,
-                           abs_tol=0.05):
+    for check in check_cfa(json_name_input, data_name_input, abs_tol=0.05):
         assert check >= THRESHOLD
 
 
-def test_13_cfa():
+def test_13_cfa():  # noqa: D103
 
-    json_name_input = 'test13'
-    data_name_input = 'test13'
+    json_name_input = "test13"
+    data_name_input = "test13"
 
-    for check in check_cfa(json_name_input,
-                           data_name_input,
-                           index_col=0,
-                           is_cov=True,
-                           n_obs=64,
-                           rel_tol=0.1):
+    for check in check_cfa(
+        json_name_input,
+        data_name_input,
+        index_col=0,
+        is_cov=True,
+        n_obs=64,
+        rel_tol=0.1,
+    ):
         assert check >= THRESHOLD
 
 
-def test_14_cfa():
+def test_14_cfa():  # noqa: D103
 
-    json_name_input = 'test14'
-    data_name_input = 'test14'
+    json_name_input = "test14"
+    data_name_input = "test14"
 
-    for check in check_cfa(json_name_input,
-                           data_name_input,
-                           index_col=0,
-                           rel_tol=0.1):
+    for check in check_cfa(json_name_input, data_name_input, index_col=0, rel_tol=0.1):
         assert check >= THRESHOLD
 
 
 @raises(ValueError)
-def test_14_cfa_no_model():
+def test_14_cfa_no_model():  # noqa: D103
 
     X = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
 
-    cfa = ConfirmatoryFactorAnalyzer('string_not_model')
+    cfa = ConfirmatoryFactorAnalyzer("string_not_model")
     cfa.fit(X)
 
 
 @raises(AssertionError)
-def test_14_cfa_bad_bounds():
+def test_14_cfa_bad_bounds():  # noqa: D103
 
     X = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
 
@@ -85,28 +83,29 @@ def test_14_cfa_bad_bounds():
 
 
 @raises(ValueError)
-def test_14_cfa_cov_with_no_obs():
+def test_14_cfa_cov_with_no_obs():  # noqa: D103
 
     ConfirmatoryFactorAnalyzer(is_cov_matrix=True)
 
 
-class TestModelSpecificationParser:
-
-    def test_model_spec_str(self):
+class TestModelSpecificationParser:  # noqa: D101
+    def test_model_spec_str(self):  # noqa: D102
 
         ms = ModelSpecification(np.array([[0, 0, 0]]), 3, 1)
-        assert str(ms).startswith('<ModelSpecification object at ')
+        assert str(ms).startswith("<ModelSpecification object at ")
 
-    def test_model_spec_as_dict(self):
+    def test_model_spec_as_dict(self):  # noqa: D102
 
         loadings = np.array([[0, 0, 0]])
         n_factors = 3
         n_variables = 1
         ms = ModelSpecification(loadings, n_factors, n_variables)
 
-        expected = {'loadings': loadings,
-                    'n_variables': n_variables,
-                    'n_factors': n_factors}
+        expected = {
+            "loadings": loadings,
+            "n_variables": n_variables,
+            "n_factors": n_factors,
+        }
         new_dict = ms.get_model_specification_as_dict()
         for key, value in expected.items():
             assert key in new_dict
@@ -115,7 +114,7 @@ def test_model_spec_as_dict(self):
             else:
                 eq_(new_dict[key], value)
 
-    def test_model_spec_parser_from_dict_none(self):
+    def test_model_spec_parser_from_dict_none(self):  # noqa: D102
         X = np.array([[0, 0, 0]])
         ms = ModelSpecificationParser.parse_model_specification_from_dict(X, None)
         assert isinstance(ms, ModelSpecification)
@@ -124,11 +123,11 @@ def test_model_spec_parser_from_dict_none(self):
         assert_array_equal(ms.loadings, np.ones((3, 3), dtype=int))
 
     @raises(ValueError)
-    def test_model_spec_parser_from_dict_error(self):
+    def test_model_spec_parser_from_dict_error(self):  # noqa: D102
         X = np.array([[0, 0, 0]])
-        ModelSpecificationParser.parse_model_specification_from_dict(X, 'not_a_model')
+        ModelSpecificationParser.parse_model_specification_from_dict(X, "not_a_model")
 
-    def test_model_spec_parser_from_array_none(self):
+    def test_model_spec_parser_from_array_none(self):  # noqa: D102
         X = np.array([[0, 0, 0]])
         ms = ModelSpecificationParser.parse_model_specification_from_array(X, None)
         assert isinstance(ms, ModelSpecification)
@@ -136,7 +135,7 @@ def test_model_spec_parser_from_array_none(self):
         eq_(ms.n_variables, 3)
         assert_array_equal(ms.loadings, np.ones((3, 3), dtype=int))
 
-    def test_model_spec_parser_from_array(self):
+    def test_model_spec_parser_from_array(self):  # noqa: D102
         X = np.array([[0, 0, 0]])
         spec = np.ones((3, 3), dtype=int)
         ms = ModelSpecificationParser.parse_model_specification_from_array(X, spec)
@@ -145,7 +144,7 @@ def test_model_spec_parser_from_array(self):
         eq_(ms.n_variables, 3)
         assert_array_equal(ms.loadings, np.ones((3, 3), dtype=int))
 
-    def test_model_spec_parser_from_frame(self):
+    def test_model_spec_parser_from_frame(self):  # noqa: D102
         X = np.array([[0, 0, 0]])
         spec = pd.DataFrame(np.ones((3, 3), dtype=int))
         ms = ModelSpecificationParser.parse_model_specification_from_array(X, spec)
@@ -155,6 +154,6 @@ def test_model_spec_parser_from_frame(self):
         assert_array_equal(ms.loadings, np.ones((3, 3), dtype=int))
 
     @raises(ValueError)
-    def test_model_spec_parser_from_array_error(self):
+    def test_model_spec_parser_from_array_error(self):  # noqa: D102
         X = np.array([[0, 0, 0]])
-        ModelSpecificationParser.parse_model_specification_from_array(X, 'not_a_model')
+        ModelSpecificationParser.parse_model_specification_from_array(X, "not_a_model")
diff --git a/tests/test_expected_factor_analyzer.py b/tests/test_expected_factor_analyzer.py
index ab0fcb9..7fe5fb7 100644
--- a/tests/test_expected_factor_analyzer.py
+++ b/tests/test_expected_factor_analyzer.py
@@ -1,9 +1,10 @@
 """
-Tests for FactorAnalyzer class
+Tests for expected factor analysis.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
-:date: 10/25/2017
-:organization: ETS
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
+:author: Nitin Madnani (nmadnani@ets.org)
+:organization: Educational Testing Service
+:date: 2022-09-05
 """
 
 from factor_analyzer.test_utils import check_scenario
@@ -14,467 +15,473 @@
 THRESHOLD = 0.9
 
 
-def test_01_none_minres_3_factors_use_corr():
+def test_01_none_minres_3_factors_use_corr():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 3
-    method = 'uls'
-    rotation = 'none'
+    method = "uls"
+    rotation = "none"
 
-    for check in check_scenario(test_name, factors, method, rotation, use_corr_matrix=True):
+    for check in check_scenario(
+        test_name, factors, method, rotation, use_corr_matrix=True
+    ):
         assert check > THRESHOLD
 
 
-def test_01_none_minres_3_factors():
+def test_01_none_minres_3_factors():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 3
-    method = 'uls'
-    rotation = 'none'
+    method = "uls"
+    rotation = "none"
 
     for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_01_none_minres_2_factors():
+def test_01_none_minres_2_factors():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 2
-    method = 'uls'
-    rotation = 'none'
+    method = "uls"
+    rotation = "none"
 
-    for check in check_scenario(test_name, factors, method, rotation,
-                                check_scores=True):
+    for check in check_scenario(
+        test_name, factors, method, rotation, check_scores=True
+    ):
         assert check > THRESHOLD
 
 
-def test_01_varimax_minres_2_factors():
+def test_01_varimax_minres_2_factors():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 2
-    method = 'uls'
-    rotation = 'varimax'
+    method = "uls"
+    rotation = "varimax"
 
-    for check in check_scenario(test_name, factors, method, rotation,
-                                check_scores=True):
+    for check in check_scenario(
+        test_name, factors, method, rotation, check_scores=True
+    ):
         assert check > THRESHOLD
 
 
-def test_01_promax_minres_2_factors():
+def test_01_promax_minres_2_factors():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 2
-    method = 'uls'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name, factors, method, rotation,
-                                ignore_value=True,
-                                ignore_communalities=True,
-                                check_structure=True,
-                                check_scores=True):
+    method = "uls"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+        check_structure=True,
+        check_scores=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_01_varimax_ml_3_factors():
+def test_01_varimax_ml_3_factors():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 3
-    method = 'ml'
-    rotation = 'varimax'
+    method = "ml"
+    rotation = "varimax"
 
     for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_01_promax_ml_3_factors():
+def test_01_promax_ml_3_factors():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 3
-    method = 'ml'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "ml"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_01_promax_ml_3_factors_use_corr():
+def test_01_promax_ml_3_factors_use_corr():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 3
-    method = 'ml'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                use_corr_matrix=True,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "ml"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        use_corr_matrix=True,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_01_promax_uls_3_factors():
+def test_01_promax_uls_3_factors():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 3
-    method = 'uls'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "uls"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_02_none_minres_3_factors():
+def test_02_none_minres_3_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 3
-    method = 'uls'
-    rotation = 'none'
+    method = "uls"
+    rotation = "none"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_02_varimax_minres_2_factors():
+def test_02_varimax_minres_2_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 2
-    method = 'uls'
-    rotation = 'varimax'
+    method = "uls"
+    rotation = "varimax"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                check_scores=True):
+    for check in check_scenario(
+        test_name, factors, method, rotation, check_scores=True
+    ):
         assert check > THRESHOLD
 
 
-def test_02_promax_minres_2_factors():
+def test_02_promax_minres_2_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 2
-    method = 'uls'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name, factors, method, rotation,
-                                ignore_value=True,
-                                ignore_communalities=True,
-                                check_structure=True,
-                                check_scores=True):
+    method = "uls"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+        check_structure=True,
+        check_scores=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_03_none_minres_3_factors():
+def test_03_none_minres_3_factors():  # noqa: D103
 
-    test_name = 'test03'
+    test_name = "test03"
     factors = 3
-    method = 'uls'
-    rotation = 'none'
+    method = "uls"
+    rotation = "none"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_03_promax_minres_3_factors():
+def test_03_promax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test03'
+    test_name = "test03"
     factors = 3
-    method = 'uls'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "uls"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_03_none_ml_2_factors():
+def test_03_none_ml_2_factors():  # noqa: D103
 
-    test_name = 'test03'
+    test_name = "test03"
     factors = 2
-    method = 'ml'
-    rotation = 'none'
+    method = "ml"
+    rotation = "none"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_04_varimax_ml_2_factors():
+def test_04_varimax_ml_2_factors():  # noqa: D103
 
-    test_name = 'test04'
+    test_name = "test04"
     factors = 2
-    method = 'ml'
-    rotation = 'varimax'
+    method = "ml"
+    rotation = "varimax"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_04_varimax_minres_3_factors():
+def test_04_varimax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test04'
+    test_name = "test04"
     factors = 3
-    method = 'uls'
-    rotation = 'varimax'
+    method = "uls"
+    rotation = "varimax"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_05_varimax_minres_3_factors():
+def test_05_varimax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test05'
+    test_name = "test05"
     factors = 3
-    method = 'uls'
-    rotation = 'varimax'
+    method = "uls"
+    rotation = "varimax"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_05_none_ml_2_factors():
+def test_05_none_ml_2_factors():  # noqa: D103
 
-    test_name = 'test05'
+    test_name = "test05"
     factors = 2
-    method = 'ml'
-    rotation = 'none'
+    method = "ml"
+    rotation = "none"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_06_promax_ml_3_factors():
+def test_06_promax_ml_3_factors():  # noqa: D103
 
-    test_name = 'test06'
+    test_name = "test06"
     factors = 3
-    method = 'ml'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "ml"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_06_none_minres_2_factors():
+def test_06_none_minres_2_factors():  # noqa: D103
 
-    test_name = 'test06'
+    test_name = "test06"
     factors = 2
-    method = 'uls'
-    rotation = 'none'
+    method = "uls"
+    rotation = "none"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_07_varimax_minres_3_factors():
+def test_07_varimax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test07'
+    test_name = "test07"
     factors = 3
-    method = 'uls'
-    rotation = 'varimax'
+    method = "uls"
+    rotation = "varimax"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_07_varimax_ml_3_factors():
+def test_07_varimax_ml_3_factors():  # noqa: D103
 
-    test_name = 'test07'
+    test_name = "test07"
     factors = 3
-    method = 'ml'
-    rotation = 'varimax'
+    method = "ml"
+    rotation = "varimax"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_08_promax_ml_3_factors():
+def test_08_promax_ml_3_factors():  # noqa: D103
 
-    test_name = 'test08'
+    test_name = "test08"
     factors = 3
-    method = 'ml'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "ml"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_08_none_ml_2_factors():
+def test_08_none_ml_2_factors():  # noqa: D103
 
-    test_name = 'test08'
+    test_name = "test08"
     factors = 2
-    method = 'ml'
-    rotation = 'none'
+    method = "ml"
+    rotation = "none"
 
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation):
+    for check in check_scenario(test_name, factors, method, rotation):
         assert check > THRESHOLD
 
 
-def test_09_promax_ml_3_factors():
+def test_09_promax_ml_3_factors():  # noqa: D103
 
-    test_name = 'test09'
+    test_name = "test09"
     factors = 3
-    method = 'ml'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "ml"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_09_promax_minres_2_factors():
+def test_09_promax_minres_2_factors():  # noqa: D103
 
-    test_name = 'test09'
+    test_name = "test09"
     factors = 2
-    method = 'uls'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "uls"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_10_none_ml_3_factors():
+def test_10_none_ml_3_factors():  # noqa: D103
 
-    test_name = 'test10'
+    test_name = "test10"
     factors = 3
-    method = 'ml'
-    rotation = 'none'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "ml"
+    rotation = "none"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_10_varimax_minres_3_factors():
+def test_10_varimax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test10'
+    test_name = "test10"
     factors = 3
-    method = 'uls'
-    rotation = 'varimax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "uls"
+    rotation = "varimax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_10_promax_minres_3_factors():
+def test_10_promax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test10'
+    test_name = "test10"
     factors = 3
-    method = 'uls'
-    rotation = 'promax'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "uls"
+    rotation = "promax"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
 
-def test_02_none_principal():
+def test_02_none_principal():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 3
-    method = 'principal'
-    rotation = 'none'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                ignore_value=True,
-                                ignore_communalities=True):
+    method = "principal"
+    rotation = "none"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
         assert check > THRESHOLD
 
-def test_15_none_principal():
 
-    test_name = 'test15'
+def test_15_none_principal():  # noqa: D103
+
+    test_name = "test15"
     factors = 20
-    method = 'principal'
-    rotation = 'none'
-    svd_method = 'lapack'
-
-    for check in check_scenario(test_name,
-                                factors,
-                                method,
-                                rotation,
-                                svd_method=svd_method,
-                                ignore_value=True,
-                                ignore_communalities=True):
-        assert check > THRESHOLD
\ No newline at end of file
+    method = "principal"
+    rotation = "none"
+    svd_method = "lapack"
+
+    for check in check_scenario(
+        test_name,
+        factors,
+        method,
+        rotation,
+        svd_method=svd_method,
+        ignore_value=True,
+        ignore_communalities=True,
+    ):
+        assert check > THRESHOLD
diff --git a/tests/test_expected_rotator.py b/tests/test_expected_rotator.py
index 243c21f..f7d4c58 100644
--- a/tests/test_expected_rotator.py
+++ b/tests/test_expected_rotator.py
@@ -1,9 +1,9 @@
-"""
-Tests for Rotator class
+"""Tests for ``Rotator`` class.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
-:date: 10/25/2017
-:organization: ETS
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
+:author: Nitin Madnani (nmadnani@ets.org)
+:organization: Educational Testing Service
+:date: 2022-09-05
 """
 
 from factor_analyzer.test_utils import check_rotation
@@ -12,240 +12,244 @@
 THRESHOLD = 0.95
 
 
-def test_01_varimax_minres_2_factors():
+def test_01_varimax_minres_2_factors():  # noqa: D103
 
-    test_name = 'test01'
+    test_name = "test01"
     factors = 2
-    method = 'uls'
-    rotation = 'varimax'
+    method = "uls"
+    rotation = "varimax"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_02_varimax_minres_3_factors():
+def test_02_varimax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 3
-    method = 'uls'
-    rotation = 'varimax'
+    method = "uls"
+    rotation = "varimax"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_02_oblimax_minres_3_factors():
+def test_02_oblimax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 3
-    method = 'uls'
-    rotation = 'oblimax'
+    method = "uls"
+    rotation = "oblimax"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_02_oblimin_minres_3_factors():
+def test_02_oblimin_minres_3_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 3
-    method = 'uls'
-    rotation = 'oblimin'
+    method = "uls"
+    rotation = "oblimin"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_02_quartimax_minres_3_factors():
+def test_02_quartimax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 3
-    method = 'uls'
-    rotation = 'quartimax'
+    method = "uls"
+    rotation = "quartimax"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_02_quartimin_minres_3_factors():
+def test_02_quartimin_minres_3_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 3
-    method = 'uls'
-    rotation = 'quartimin'
+    method = "uls"
+    rotation = "quartimin"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_04_oblimax_minres_3_factors():
+def test_04_oblimax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test04'
+    test_name = "test04"
     factors = 3
-    method = 'uls'
-    rotation = 'oblimax'
+    method = "uls"
+    rotation = "oblimax"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_04_oblimin_minres_3_factors():
+def test_04_oblimin_minres_3_factors():  # noqa: D103
 
-    test_name = 'test04'
+    test_name = "test04"
     factors = 3
-    method = 'uls'
-    rotation = 'oblimin'
+    method = "uls"
+    rotation = "oblimin"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_04_quartimax_minres_3_factors():
+def test_04_quartimax_minres_3_factors():  # noqa: D103
 
-    test_name = 'test04'
+    test_name = "test04"
     factors = 3
-    method = 'uls'
-    rotation = 'quartimax'
+    method = "uls"
+    rotation = "quartimax"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_04_quartimin_minres_3_factors():
+def test_04_quartimin_minres_3_factors():  # noqa: D103
 
-    test_name = 'test04'
+    test_name = "test04"
     factors = 3
-    method = 'uls'
-    rotation = 'quartimin'
+    method = "uls"
+    rotation = "quartimin"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_07_oblimax_minres_2_factors():
+def test_07_oblimax_minres_2_factors():  # noqa: D103
 
-    test_name = 'test07'
+    test_name = "test07"
     factors = 2
-    method = 'uls'
-    rotation = 'oblimax'
+    method = "uls"
+    rotation = "oblimax"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_07_oblimin_minres_3_factors():
+def test_07_oblimin_minres_3_factors():  # noqa: D103
 
-    test_name = 'test07'
+    test_name = "test07"
     factors = 2
-    method = 'uls'
-    rotation = 'oblimin'
+    method = "uls"
+    rotation = "oblimin"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_07_quartimax_minres_2_factors():
+def test_07_quartimax_minres_2_factors():  # noqa: D103
 
-    test_name = 'test07'
+    test_name = "test07"
     factors = 2
-    method = 'uls'
-    rotation = 'quartimax'
+    method = "uls"
+    rotation = "quartimax"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_07_quartimin_minres_2_factors():
+def test_07_quartimin_minres_2_factors():  # noqa: D103
 
-    test_name = 'test07'
+    test_name = "test07"
     factors = 2
-    method = 'uls'
-    rotation = 'quartimin'
+    method = "uls"
+    rotation = "quartimin"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_07_equamax_minres_2_factors():
+def test_07_equamax_minres_2_factors():  # noqa: D103
 
-    test_name = 'test07'
+    test_name = "test07"
     factors = 2
-    method = 'uls'
-    rotation = 'equamax'
+    method = "uls"
+    rotation = "equamax"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check > THRESHOLD
 
 
-def test_07_oblimin_minres_2_factors_gamma():
+def test_07_oblimin_minres_2_factors_gamma():  # noqa: D103
 
-    test_name = 'test07_gamma'
+    test_name = "test07_gamma"
     factors = 2
-    method = 'uls'
-    rotation = 'oblimin'
+    method = "uls"
+    rotation = "oblimin"
     gamma = 0.5
 
     check = check_rotation(test_name, factors, method, rotation, gamma=gamma)
     assert check > THRESHOLD
 
 
-def test_02_geomin_obl_minres_2_factors():
+def test_02_geomin_obl_minres_2_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 2
-    method = 'uls'
-    rotation = 'geomin_obl'
+    method = "uls"
+    rotation = "geomin_obl"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check == 1
 
-def test_02_geomin_obl_minres_3_factors():
 
-    test_name = 'test02'
+def test_02_geomin_obl_minres_3_factors():  # noqa: D103
+
+    test_name = "test02"
     factors = 3
-    method = 'uls'
-    rotation = 'geomin_obl'
+    method = "uls"
+    rotation = "geomin_obl"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check == 1
 
-def test_02_geomin_obl_ml_3_factors():
 
-    test_name = 'test02'
+def test_02_geomin_obl_ml_3_factors():  # noqa: D103
+
+    test_name = "test02"
     factors = 3
-    method = 'ml'
-    rotation = 'geomin_obl'
+    method = "ml"
+    rotation = "geomin_obl"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check == 1
 
-def test_02_geomin_ort_minres_2_factors():
 
-    test_name = 'test02'
+def test_02_geomin_ort_minres_2_factors():  # noqa: D103
+
+    test_name = "test02"
     factors = 2
-    method = 'uls'
-    rotation = 'geomin_ort'
+    method = "uls"
+    rotation = "geomin_ort"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check == 1
 
-def test_02_geomin_ort_minres_3_factors():
 
-    test_name = 'test02'
+def test_02_geomin_ort_minres_3_factors():  # noqa: D103
+
+    test_name = "test02"
     factors = 3
-    method = 'uls'
-    rotation = 'geomin_ort'
+    method = "uls"
+    rotation = "geomin_ort"
 
     check = check_rotation(test_name, factors, method, rotation)
     assert check == 1
 
 
-def test_02_geomin_ort_ml_3_factors():
+def test_02_geomin_ort_ml_3_factors():  # noqa: D103
 
-    test_name = 'test02'
+    test_name = "test02"
     factors = 3
-    method = 'ml'
-    rotation = 'geomin_ort'
+    method = "ml"
+    rotation = "geomin_ort"
 
     check = check_rotation(test_name, factors, method, rotation)
-    assert check == 1
\ No newline at end of file
+    assert check == 1
diff --git a/tests/test_factor_analyzer.py b/tests/test_factor_analyzer.py
index 9302468..24cda70 100644
--- a/tests/test_factor_analyzer.py
+++ b/tests/test_factor_analyzer.py
@@ -1,41 +1,41 @@
 """
-Tests for FactorAnalyzer class
+Tests for the ``FactorAnalyzer`` class.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
-:date: 10/25/2017
-:organization: ETS
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
+:author: Nitin Madnani (nmadnani@ets.org)
+:organization: Educational Testing Service
+:date: 2022-09-05
 """
 
 import numpy as np
 import pandas as pd
-
-from nose.tools import raises
+from nose.tools import assert_almost_equal, raises
 from numpy.testing import assert_array_almost_equal
-from pandas.util.testing import assert_almost_equal
-
 from sklearn.model_selection import GridSearchCV
 from sklearn.pipeline import make_pipeline
 from sklearn.tree import DecisionTreeClassifier
 
+from factor_analyzer.factor_analyzer import (
+    FactorAnalyzer,
+    calculate_bartlett_sphericity,
+    calculate_kmo,
+)
 from factor_analyzer.utils import smc
-from factor_analyzer.factor_analyzer import FactorAnalyzer
-from factor_analyzer.factor_analyzer import (calculate_kmo,
-                                             calculate_bartlett_sphericity)
 
 
-def test_calculate_bartlett_sphericity():
+def test_calculate_bartlett_sphericity():  # noqa: D103
 
-    path = 'tests/data/test01.csv'
+    path = "tests/data/test01.csv"
     data = pd.read_csv(path)
     s, p = calculate_bartlett_sphericity(data.values)
 
-    assert_almost_equal(s, 14185.00286)
-    assert_almost_equal(p, 0)
+    assert_almost_equal(s, 14185, places=2)
+    assert_almost_equal(p, 0, places=2)
 
 
-def test_calculate_kmo():
+def test_calculate_kmo():  # noqa: D103
 
-    path = 'tests/data/test02.csv'
+    path = "tests/data/test02.csv"
     data = pd.read_csv(path)
 
     expected_overall = 0.81498
@@ -45,139 +45,167 @@ def test_calculate_kmo():
         0.56004925345997,
         0.7000330131087749,
         0.7054455920793854,
-        0.829063299715461, 
+        0.829063299715461,
         0.8484249727243623,
         0.8635016393452357,
         0.8411432114912387,
         0.8770763964694772,
-        0.8392720039048369]
+        0.8392720039048369,
+    ]
 
     expected_by_item = np.array(values)
 
-    (kmo_by_item,
-     kmo_overall) = calculate_kmo(data.values)
+    (kmo_by_item, kmo_overall) = calculate_kmo(data.values)
 
-    assert_almost_equal(kmo_by_item, expected_by_item)
-    assert_almost_equal(kmo_overall, expected_overall)
+    assert_array_almost_equal(kmo_by_item, expected_by_item)
+    assert_almost_equal(kmo_overall, expected_overall, places=5)
 
 
-def test_gridsearch():
+def test_gridsearch():  # noqa: D103
     # make sure this doesn't fail
 
     X = pd.DataFrame(np.random.randn(1000).reshape(100, 10))
     y = pd.Series(np.random.choice([1, 0], size=100))
 
-    grid = {'factoranalyzer__n_factors': [5, 7],
-            'factoranalyzer__rotation': [None, 'varimax'],
-            'decisiontreeclassifier__max_depth': [2, 5]}
+    grid = {
+        "factoranalyzer__n_factors": [5, 7],
+        "factoranalyzer__rotation": [None, "varimax"],
+        "decisiontreeclassifier__max_depth": [2, 5],
+    }
 
     fa = FactorAnalyzer()
     decisiontree = DecisionTreeClassifier(random_state=123)
     pipe = make_pipeline(fa, decisiontree)
 
-    gridsearch = GridSearchCV(pipe,
-                              grid,
-                              scoring='f1',
-                              cv=3,
-                              verbose=0)
+    gridsearch = GridSearchCV(pipe, grid, scoring="f1", cv=3, verbose=0)
     gridsearch.fit(X, y)
 
 
-class TestFactorAnalyzer:
-
-    def test_analyze_weights(self):
+class TestFactorAnalyzer:  # noqa: D101
+    def test_analyze_weights(self):  # noqa: D102
 
-        data = pd.DataFrame({'A': [2, 4, 5, 6, 8, 9],
-                             'B': [4, 8, 9, 10, 16, 18],
-                             'C': [6, 12, 15, 12, 26, 27]})
+        data = pd.DataFrame(
+            {
+                "A": [2, 4, 5, 6, 8, 9],
+                "B": [4, 8, 9, 10, 16, 18],
+                "C": [6, 12, 15, 12, 26, 27],
+            }
+        )
 
         fa = FactorAnalyzer(rotation=None)
         fa.fit(data)
         _ = fa.transform(data)
-        expected_weights = np.array(([[0.335363, -2.725096,  0.],
-                                      [0.339166, -0.293888,  0.],
-                                      [0.334446,  3.030608, -0.]]))
+        expected_weights = np.array(
+            (
+                [
+                    [0.335363, -2.725096, 0.0],
+                    [0.339166, -0.293888, 0.0],
+                    [0.334446, 3.030608, -0.0],
+                ]
+            )
+        )
         assert_array_almost_equal(expected_weights, fa.weights_)
 
-    def test_analyze_impute_mean(self):
+    def test_analyze_impute_mean(self):  # noqa: D102
 
-        data = pd.DataFrame({'A': [2, 4, 5, 6, 8, 9],
-                             'B': [4, 8, np.nan, 10, 16, 18],
-                             'C': [6, 12, 15, 12, 26, 27]})
+        data = pd.DataFrame(
+            {
+                "A": [2, 4, 5, 6, 8, 9],
+                "B": [4, 8, np.nan, 10, 16, 18],
+                "C": [6, 12, 15, 12, 26, 27],
+            }
+        )
 
         expected = data.copy()
         expected.iloc[2, 1] = np.mean([4, 8, 10, 16, 18])
         expected_corr = expected.corr()
         expected_corr = expected_corr.values
 
-        fa = FactorAnalyzer(rotation=None, impute='mean', n_factors=1)
+        fa = FactorAnalyzer(rotation=None, impute="mean", n_factors=1)
         fa.fit(data)
         assert_array_almost_equal(fa.corr_, expected_corr)
 
-    def test_analyze_impute_median(self):
+    def test_analyze_impute_median(self):  # noqa: D102
 
-        data = pd.DataFrame({'A': [2, 4, 5, 6, 8, 9],
-                             'B': [4, 8, np.nan, 10, 16, 18],
-                             'C': [6, 12, 15, 12, 26, 27]})
+        data = pd.DataFrame(
+            {
+                "A": [2, 4, 5, 6, 8, 9],
+                "B": [4, 8, np.nan, 10, 16, 18],
+                "C": [6, 12, 15, 12, 26, 27],
+            }
+        )
 
         expected_corr = pd.DataFrame(
             [
-                [1.      , 0.987534, 0.946014],
-                [0.987534, 1.      , 0.983488],
-                [0.946014, 0.983488, 1.      ]
-        ]).values
+                [1.0, 0.987534, 0.946014],
+                [0.987534, 1.0, 0.983488],
+                [0.946014, 0.983488, 1.0],
+            ]
+        ).values
 
-        fa = FactorAnalyzer(rotation=None, impute='median', n_factors=1)
+        fa = FactorAnalyzer(rotation=None, impute="median", n_factors=1)
         fa.fit(data)
         assert_array_almost_equal(expected_corr, fa.corr_)
 
-    def test_analyze_impute_drop(self):
+    def test_analyze_impute_drop(self):  # noqa: D102
 
-        data = pd.DataFrame({'A': [2, 4, 5, 6, 8, 9],
-                             'B': [4, 8, np.nan, 10, 16, 18],
-                             'C': [6, 12, 15, 12, 26, 27]})
+        data = pd.DataFrame(
+            {
+                "A": [2, 4, 5, 6, 8, 9],
+                "B": [4, 8, np.nan, 10, 16, 18],
+                "C": [6, 12, 15, 12, 26, 27],
+            }
+        )
 
         expected = data.copy()
         expected = expected.dropna()
         expected_corr = expected.corr()
         expected_corr = expected_corr.values
 
-        fa = FactorAnalyzer(rotation=None, impute='drop', n_factors=1)
+        fa = FactorAnalyzer(rotation=None, impute="drop", n_factors=1)
         fa.fit(data)
         assert_array_almost_equal(fa.corr_, expected_corr)
 
     @raises(ValueError)
-    def test_analyze_bad_svd_method(self):
-        fa = FactorAnalyzer(svd_method='foo')
+    def test_analyze_bad_svd_method(self):  # noqa: D102
+        fa = FactorAnalyzer(svd_method="foo")
         fa.fit(np.random.randn(500).reshape(100, 5))
 
     @raises(ValueError)
-    def test_analyze_impute_value_error(self):
-        fa = FactorAnalyzer(rotation=None, impute='blah', n_factors=1)
+    def test_analyze_impute_value_error(self):  # noqa: D102
+        fa = FactorAnalyzer(rotation=None, impute="blah", n_factors=1)
         fa.fit(np.random.randn(500).reshape(100, 5))
 
     @raises(ValueError)
-    def test_analyze_rotation_value_error(self):
-        fa = FactorAnalyzer(rotation='blah', n_factors=1)
+    def test_analyze_rotation_value_error(self):  # noqa: D102
+        fa = FactorAnalyzer(rotation="blah", n_factors=1)
         fa.fit(np.random.randn(500).reshape(100, 5))
 
     @raises(ValueError)
-    def test_analyze_infinite(self):
-
-        data = pd.DataFrame({'A': [1.0, 0.4, 0.5],
-                             'B': [0.4, 1.0, float('inf')],
-                             'C': [0.5, float('inf'), 1.0]},
-                            index=['A', 'B', 'C'])
-
-        fa = FactorAnalyzer(impute='drop', n_factors=1, is_corr_matrix=True)
+    def test_analyze_infinite(self):  # noqa: D102
+
+        data = pd.DataFrame(
+            {
+                "A": [1.0, 0.4, 0.5],
+                "B": [0.4, 1.0, float("inf")],
+                "C": [0.5, float("inf"), 1.0],
+            },
+            index=["A", "B", "C"],
+        )
+
+        fa = FactorAnalyzer(impute="drop", n_factors=1, is_corr_matrix=True)
         fa.fit(data)
 
-    def test_smc_is_r_squared(self):
+    def test_smc_is_r_squared(self):  # noqa: D102
         # test that SMC is roughly equivalent to R-squared values.
 
-        data = pd.DataFrame({'A': [10.5, 20.1, 30.2, 40.1, 50.3],
-                             'B': [62, 71, 83, 91, 15],
-                             'C': [0.45, 0.90, 0.22, 0.34, .045]})
+        data = pd.DataFrame(
+            {
+                "A": [10.5, 20.1, 30.2, 40.1, 50.3],
+                "B": [62, 71, 83, 91, 15],
+                "C": [0.45, 0.90, 0.22, 0.34, 0.045],
+            }
+        )
 
         expected_r2 = [0.478330, 0.196223, 0.484519]
         expected_r2 = np.array(expected_r2)
@@ -186,9 +214,9 @@ def test_smc_is_r_squared(self):
 
         assert_array_almost_equal(smc_result, expected_r2)
 
-    def test_factor_variance(self):
+    def test_factor_variance(self):  # noqa: D102
 
-        path = 'tests/data/test01.csv'
+        path = "tests/data/test01.csv"
         data = pd.read_csv(path)
 
         fa = FactorAnalyzer(n_factors=3, rotation=None)
@@ -198,11 +226,11 @@ def test_factor_variance(self):
         n_rows = loadings.shape[0]
 
         # calculate variance
-        loadings = loadings ** 2
+        loadings = loadings**2
         variance = np.sum(loadings, axis=0)
 
         # calculate proportional variance
         proportional_variance_expected = variance / n_rows
         proportional_variance = fa.get_factor_variance()[1]
 
-        assert_almost_equal(proportional_variance_expected, proportional_variance)
+        assert_array_almost_equal(proportional_variance_expected, proportional_variance)
diff --git a/tests/test_utils.py b/tests/test_utils.py
index e389687..ccc0e26 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -1,32 +1,35 @@
 """
 Tests utilities.
 
-:author: Jeremy Biggs (jbiggs@ets.org)
-:date: 02/15/2019
-:organization: ETS
+:author: Jeremy Biggs (jeremy.m.biggs@gmail.com)
+:author: Nitin Madnani (nmadnani@ets.org)
+:organization: Educational Testing Service
+:date: 2022-09-05
 """
 
 import numpy as np
 import pandas as pd
-from factor_analyzer.utils import (covariance_to_correlation,
-                                   duplication_matrix,
-                                   duplication_matrix_pre_post,
-                                   commutation_matrix,
-                                   fill_lower_diag,
-                                   get_first_idxs_from_values,
-                                   get_free_parameter_idxs,
-                                   get_symmetric_lower_idxs,
-                                   get_symmetric_upper_idxs,
-                                   merge_variance_covariance,
-                                   partial_correlations,
-                                   unique_elements)
-
 from nose.tools import eq_, raises
 from numpy.testing import assert_array_equal
 from pandas.util.testing import assert_almost_equal
 
+from factor_analyzer.utils import (
+    commutation_matrix,
+    covariance_to_correlation,
+    duplication_matrix,
+    duplication_matrix_pre_post,
+    fill_lower_diag,
+    get_first_idxs_from_values,
+    get_free_parameter_idxs,
+    get_symmetric_lower_idxs,
+    get_symmetric_upper_idxs,
+    merge_variance_covariance,
+    partial_correlations,
+    unique_elements,
+)
+
 
-def test_unique_elements():
+def test_unique_elements():  # noqa: D103
 
     expected = [1, 2, 3, 4, 5]
 
@@ -35,16 +38,14 @@ def test_unique_elements():
     eq_(output, expected)
 
 
-def test_fill_lower_diag():
+def test_fill_lower_diag():  # noqa: D103
 
-    expected = np.array([[0, 0, 0],
-                         [1, 0, 0],
-                         [2, 3, 0]])
+    expected = np.array([[0, 0, 0], [1, 0, 0], [2, 3, 0]])
     output = fill_lower_diag([1, 2, 3])
     assert_array_equal(output, expected)
 
 
-def test_merge_variance_covariance_no_covariance():
+def test_merge_variance_covariance_no_covariance():  # noqa: D103
 
     expected = np.eye(4)
 
@@ -54,36 +55,29 @@ def test_merge_variance_covariance_no_covariance():
     assert_array_equal(output, expected)
 
 
-def test_merge_variance_covariance():
+def test_merge_variance_covariance():  # noqa: D103
 
-    expected = [[1, .25, .45],
-                [.25, 1, .35],
-                [.45, .35, 1]]
+    expected = [[1, 0.25, 0.45], [0.25, 1, 0.35], [0.45, 0.35, 1]]
     expected = np.array(expected)
 
     x = np.array([1, 1, 1])
-    y = np.array([.25, .45, .35])
+    y = np.array([0.25, 0.45, 0.35])
 
     output = merge_variance_covariance(x, y)
     assert_array_equal(output, expected)
 
 
-def test_get_first_idxs_from_values():
+def test_get_first_idxs_from_values():  # noqa: D103
     expected = [2, 0, 1], [0, 1, 2]
-    x = np.array([[np.nan, 1, np.nan],
-                  [np.nan, 1, 1],
-                  [1, 2, np.nan],
-                  [1, 1, 1]])
+    x = np.array([[np.nan, 1, np.nan], [np.nan, 1, 1], [1, 2, np.nan], [1, 1, 1]])
     output = get_first_idxs_from_values(x)
 
     eq_(output, expected)
 
 
-def test_get_free_parameter_idxs():
+def test_get_free_parameter_idxs():  # noqa: D103
 
-    x = np.array([[np.nan, np.nan, 1],
-                  [1, np.nan, 1],
-                  [1, 1, np.nan]])
+    x = np.array([[np.nan, np.nan, 1], [1, np.nan, 1], [1, 1, np.nan]])
 
     expected = np.array([0, 3, 4, 8])
 
@@ -91,33 +85,32 @@ def test_get_free_parameter_idxs():
     assert_array_equal(output, expected)
 
 
-def test_duplication_matrix():
+def test_duplication_matrix():  # noqa: D103
 
     output = duplication_matrix(3)
 
-    expected = np.array([[1., 0., 0., 0., 0., 0.],
-                         [0., 1., 0., 0., 0., 0.],
-                         [0., 0., 1., 0., 0., 0.],
-                         [0., 1., 0., 0., 0., 0.],
-                         [0., 0., 0., 1., 0., 0.],
-                         [0., 0., 0., 0., 1., 0.],
-                         [0., 0., 1., 0., 0., 0.],
-                         [0., 0., 0., 0., 1., 0.],
-                         [0., 0., 0., 0., 0., 1.]])
+    expected = np.array(
+        [
+            [1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
+            [0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
+            [0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
+            [0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
+            [0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
+        ]
+    )
 
     assert_array_equal(output, expected)
 
 
-def test_duplication_matrix_pre_post():
+def test_duplication_matrix_pre_post():  # noqa: D103
 
-    x = np.array([[1, 2, 3, 4],
-                  [5, 6, 7, 8],
-                  [9, 10, 11, 12],
-                  [13, 14, 15, 16]])
+    x = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
 
-    expected = np.array([[1, 5, 4],
-                         [14, 34, 20],
-                         [13, 29, 16]])
+    expected = np.array([[1, 5, 4], [14, 34, 20], [13, 29, 16]])
     expected = pd.DataFrame(expected)
 
     output = duplication_matrix_pre_post(x)
@@ -125,61 +118,61 @@ def test_duplication_matrix_pre_post():
     assert_array_equal(output, expected.values)
 
 
-def test_commutation_matrix():
+def test_commutation_matrix():  # noqa: D103
 
-    expected = np.array([[1., 0., 0., 0.],
-                         [0., 0., 1., 0.],
-                         [0., 1., 0., 0.],
-                         [0., 0., 0., 1.]])
+    expected = np.array(
+        [
+            [1.0, 0.0, 0.0, 0.0],
+            [0.0, 0.0, 1.0, 0.0],
+            [0.0, 1.0, 0.0, 0.0],
+            [0.0, 0.0, 0.0, 1.0],
+        ]
+    )
 
     output = commutation_matrix(2, 2)
 
     assert_array_equal(output, expected)
 
 
-def test_get_symmetric_lower_idxs():
+def test_get_symmetric_lower_idxs():  # noqa: D103
 
-    expected = np.array([0, 1, 2, 3, 4, 6, 7, 8,
-                         9, 12, 13, 14, 18, 19, 24])
+    expected = np.array([0, 1, 2, 3, 4, 6, 7, 8, 9, 12, 13, 14, 18, 19, 24])
 
     output = get_symmetric_lower_idxs(5)
 
     assert_array_equal(output, expected)
 
 
-def test_get_symmetric_lower_idxs_no_diag():
+def test_get_symmetric_lower_idxs_no_diag():  # noqa: D103
 
-    expected = np.array([5, 10, 11, 15, 16,
-                         17, 20, 21, 22, 23])
+    expected = np.array([5, 10, 11, 15, 16, 17, 20, 21, 22, 23])
 
     output = get_symmetric_lower_idxs(5, diag=False)
 
     assert_array_equal(output, expected)
 
 
-def test_get_symmetric_upper_idxs():
+def test_get_symmetric_upper_idxs():  # noqa: D103
 
-    expected = np.array([0, 5, 10, 15, 20, 6, 11, 16,
-                         21, 12, 17, 22, 18, 23, 24])
+    expected = np.array([0, 5, 10, 15, 20, 6, 11, 16, 21, 12, 17, 22, 18, 23, 24])
 
     output = get_symmetric_upper_idxs(5)
 
     assert_array_equal(output, expected)
 
 
-def test_get_symmetric_upper_idxs_no_diag():
+def test_get_symmetric_upper_idxs_no_diag():  # noqa: D103
 
-    expected = np.array([1, 2, 7, 3, 8, 13,
-                         4, 9, 14, 19])
+    expected = np.array([1, 2, 7, 3, 8, 13, 4, 9, 14, 19])
 
     output = get_symmetric_upper_idxs(5, diag=False)
 
     assert_array_equal(output, expected)
 
 
-def test_covariance_to_correlation():
+def test_covariance_to_correlation():  # noqa: D103
 
-    path = 'tests/data/test02.csv'
+    path = "tests/data/test02.csv"
     data = pd.read_csv(path)
 
     expected_corr = data.corr().values
@@ -190,65 +183,61 @@ def test_covariance_to_correlation():
 
 
 @raises(ValueError)
-def test_covariance_to_correlation_value_error():
+def test_covariance_to_correlation_value_error():  # noqa: D103
 
-    covariance_to_correlation(np.array([[23, 12, 23],
-                                        [42, 25, 21]]))
+    covariance_to_correlation(np.array([[23, 12, 23], [42, 25, 21]]))
 
 
-def test_partial_correlations():
+def test_partial_correlations():  # noqa: D103
 
-    data = pd.DataFrame([[12, 14, 15],
-                         [24, 12, 52],
-                         [35, 12, 41],
-                         [23, 12, 42]])
+    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]]
+    expected = [
+        [1.0, -0.730955, -0.50616],
+        [-0.730955, 1.0, -0.928701],
+        [-0.50616, -0.928701, 1.0],
+    ]
 
-    expected = pd.DataFrame(expected,
-                            columns=[0, 1, 2],
-                            index=[0, 1, 2])
+    expected = pd.DataFrame(expected, columns=[0, 1, 2], index=[0, 1, 2])
 
     result = partial_correlations(data)
     assert_almost_equal(result, expected.values)
 
 
-def test_partial_correlations_num_columns_greater():
+def test_partial_correlations_num_columns_greater():  # noqa: D103
 
     # columns greater than rows
-    data = pd.DataFrame([[23, 12, 23],
-                         [42, 25, 21]])
+    data = pd.DataFrame([[23, 12, 23], [42, 25, 21]])
 
     empty_array = np.empty((3, 3))
     empty_array[:] = np.nan
     np.fill_diagonal(empty_array, 1.0)
 
-    expected = pd.DataFrame(empty_array,
-                            columns=[0, 1, 2],
-                            index=[0, 1, 2])
+    expected = pd.DataFrame(empty_array, columns=[0, 1, 2], index=[0, 1, 2])
 
     result = partial_correlations(data)
     assert_almost_equal(result, expected.values)
 
 
-def test_partial_correlations_with_zero_det():
+def test_partial_correlations_with_zero_det():  # noqa: D103
 
     # Covariance matrix that will be singular
-    data = pd.DataFrame([[10, 10, 10, 10],
-                         [12, 12, 12, 12],
-                         [15, 15, 15, 15],
-                         [20, 20, 20, 20],
-                         [11, 11, 11, 11]])
-
-    expected = [[1.0,
-                 -0.9999999999999998,
-                 -0.9999999999999998,
-                 -0.9999999999999998],
-                [-1.0000000000000004, 1.0, -1.0, -1.0],
-                [-1.0000000000000004, -1.0, 1.0, -1.0],
-                [-1.0000000000000004, -1.0, -1.0, 1.0]]
+    data = pd.DataFrame(
+        [
+            [10, 10, 10, 10],
+            [12, 12, 12, 12],
+            [15, 15, 15, 15],
+            [20, 20, 20, 20],
+            [11, 11, 11, 11],
+        ]
+    )
+
+    expected = [
+        [1.0, -0.9999999999999998, -0.9999999999999998, -0.9999999999999998],
+        [-1.0000000000000004, 1.0, -1.0, -1.0],
+        [-1.0000000000000004, -1.0, 1.0, -1.0],
+        [-1.0000000000000004, -1.0, -1.0, 1.0],
+    ]
     expected = pd.DataFrame(expected)
 
     result = partial_correlations(data)