本記事の目的

因子分析を学習していたので、備忘録的に記載します。
少しでも参考になれば嬉しいです。

因子分析とは

共通因子（以下の$f1$や$f2$）を見つけ出す手法ではなく、その根本は「因子負荷量（以下の$α_{ji}$）の値を見つけ出す手法」である（なぜか：そもそも共通因子については仮説を立てから分析するため、共通因子を見つけ出すことが目的ではなく、その係数（因子負荷量）を求めることが本来の目的である）
即ち、共通因子の想定（仮説）をある程度立ててからでないと、うまくいくものではなく、とりあえず手元にあるデータで実行しても思うような結果が出るものではない（例えば、アンケート分析の場合は、因子分析を想定したアンケート設計が必須である）
共通因子と独自因子、独自因子同士は、独立（相関はしていない）という仮定がある

回転とは

ざっくり言うと「因子負荷量を解釈しやすくするために行う変換」のこと
様々な回転方法があるが、有名なものは「バリマックス回転（直交回転）」、「プロマックス回転（斜交回転）」等である
回転を行わない場合、主成分分析のように第1成分（因子）に総合力のような共通因子（全て+になるような因子）が抽出され、第2成分以降に±が混在することで、因子の解釈が難しくなる。これを解決するのが「回転」であり、軸を動かすことで解釈しやすくしている
共通因子同士は相関しない、と仮定する考え方を「直交因子モデル」、共通因子同士の相関を仮定する考え方を「斜交因子モデル」という）
以下の回転イメージで、「回転」によって因子負荷量の解釈性が高まることを理解できる

使用するライブラリ

import matplotlib.pyplot as plt
import pandas as pd
pd.set_option('display.max_columns',None)
pd.set_option('display.max_rows',None)
import numpy as np
import japanize_matplotlib
import seaborn as sns
import warnings
warnings.simplefilter("ignore")
import math
import statsmodels.api as sm
from factor_analyzer import FactorAnalyzer as FA

使用するデータ

擬似的に作成したアンケートデータを使用していきます。簡単にアンケート実施の背景や目的を設定しておきます。
【背景】
あなたはとある動画配信サブスク「Gift（仮称）」のマーケティングを担当しています。直近1か月内に利用歴のあるランダムなユーザーに対し、顧客満足度アンケートを実施します。また、経験的にGiftは「動画コンテンツの充実感」や「アプリの機能性」により選ばれているサービスであると思われるため、これを仮説とします。
なお、アンケート対象に選ばれたユーザー全員から回答が得られたものとします。

【目的】
本満足度アンケートを通して、①：現状のサービス満足度を評価する、②：各満足度の背景にある因子を推察し、仮説と相違ないか検証する、③：年齢層別で傾向が異なるか比較する、④：①～③から今後も継続してもらうための改善示唆を出す。

回答はリッカート尺度で5段階とし、以下の設問としました。
Q1:コンテンツの充実（配信しているコンテンツ量に満足しているか？）
Q2:ジャンルの充実（幅広いジャンルに対応しているか？）
Q3:コンテンツの更新頻度（コンテンツの更新頻度に満足しているか？）
Q4:コンテンツの見つけやすさ（見たい作品はすぐに見つけられるか？）
Q5:字幕・吹替対応（字幕または吹替で見たいと思った作品の要望に応えられているか？）
Q6:同時視聴可能数（同時視聴数に満足しているか？）
Q7:視聴の質（画質・音質に満足しているか？）
Q8:価格（現在の価格に満足しているか？）

擬似データの作成

sample_size = 500

# 乱数シードの固定
np.random.seed(1234)

# 年齢を生成（乱数）
age = np.random.randint(18,70,sample_size)

# 乱数シードの固定
np.random.seed(1234)

# 年齢を生成（乱数）
age = np.random.randint(18,70,sample_size)

# 想定した因子のベースを生成
# コンテンツの豊富さ
fact1 = np.random.randint(1,6,sample_size)

# アプリの機能性
fact2 = np.random.randint(1,6,sample_size)

# 因子をベースに各質問の回答を作成
# コンテンツの充実
q1_error = np.random.randint(1,6,sample_size)
q1 = 0.7*fact1 + 0.3*q1_error

# ジャンルの充実
# 意図的に高年齢層のジャンル満足度が低下するように設定
q2_error = np.random.randint(1,6,sample_size)
ageweight = []
for a in age:
    if a>=50:
        val = np.random.randint(1,4,1)
        
    else:
        val = np.random.randint(4,6,1)
    ageweight.append(int(val))
ageweight = np.array(ageweight)
q2 = 0.3*q1 + 0.4*ageweight + 0.3*q2_error

# コンテンツの更新頻度
q3_error = np.random.randint(1,6,sample_size)
q3 = 0.6*fact1 + 0.4*q3_error

# コンテンツの見つけやすさ
q4_error = np.random.randint(1,6,sample_size)
q4 = 0.7*fact2 + 0.3*q4_error

# 字幕・吹替対応
q5_error = np.random.randint(1,6,sample_size)
q5 = 0.6*fact2 + 0.4*q5_error

# 同時視聴可能数
q6_error = np.random.randint(1,6,sample_size)
q6 = 0.5*fact2 + 0.5*q6_error

# 視聴の質
q7_error = np.random.randint(1,6,sample_size)
q7 = 0.5*fact2 + 0.5*q7_error

# 価格
# 意図的に若年齢層の価格満足度が低下するように設定
q8_error = np.random.randint(1,6,sample_size)
ageweight = []
for a in age:
    if a<30:
        val = np.random.randint(1,4,1)
        
    else:
        val = np.random.randint(4,6,1)
    ageweight.append(int(val))
ageweight = np.array(ageweight)
q8 = 0.3*fact1 + 0.1*fact2 + 0.3*ageweight + 0.3*q8_error 

data = pd.DataFrame({'年齢':age,
                      'コンテンツの充実':q1,
                      'ジャンルの充実':q2,
                      'コンテンツの更新頻度':q3,
                      'コンテンツの見つけやすさ':q4,
                      '字幕・吹替対応':q5,
                      '同時視聴可能数':q6,
                      '視聴の質':q7,
                      '価格':q8})

# 小数点を四捨五入して整数に
data = round(data,0)
data = data.astype(int)

# 要約統計量
data.describe()

年齢層別の平均満足度

20代以下（18~29歳）、30-40代（30~49歳）、50代以上（50~69歳）の3区分に分け、平均満足度を算出

data['年齢層'] = '20代以下'
data.loc[(data['年齢']>=30)&(data['年齢']<50),'年齢層'] = '30-40代'
data.loc[data['年齢']>=50,'年齢層'] = '50代以上'
data.groupby('年齢層').mean()

他の年齢層に比べ、50代以上は「ジャンルの充実」満足度、20代以下では「価格」満足度の平均値が低いとわかります。その他の項目に関しては、どの年齢層も同程度の満足度であることがわかります。
年齢層によって、サービス改善の方向性が変わってきそうです。

因子分析（Factor-analyzer）

因子数の仮定

事前に因子数を仮定する方法は、「ガットマン基準」や「スクリ―基準」等、いくつかあるようです。
ここでは「スクリ―基準」を採用し、スクリ―プロットを描画して、その推移がなだらかになる前までの因子数で仮定することとします。

【スクリ―プロット】
横軸に因子数、縦軸に固有値を取り、その推移を可視化したもの。
因子数が多くなるほど、当然モデル解釈は煩雑になるため、その因子数で観測値のデータが十分表現できる数を探す必要がある。簡単にまとめると、固有値がガクッと下がり、なだらかになり始める数（因子数）というのは、データの分散（固有値）が小さくなり始める数、即ち、その因子以降はあまり情報残ってないですよ、という数。そのため、固有値がガクッと下がり、なだらかになっていく直前の因子数を仮定する。

corr = data.drop(['年齢','年齢層'],axis=1).corr()

# 固有値、固有ベクトルを求める
eig_val, eig_vec = np.linalg.eig(corr)

# 使用するのは固有値のみ
eig_val = list(eig_val)
eig_val.sort(reverse=True)

# スクリープロットを描画
x_ = [i+1 for i in range(len(corr))]
plt.figure(figsize=(5,3))
plt.plot(x_,eig_val, marker='o')
plt.xlabel('因子数')
plt.ylabel('固有値')
plt.show()

上記のスクリ―プロットから、固有値がガクッと下がり、なだらかになり始める因子数は3であることがわかる。そのため、因子数は直前の2を仮定する（そもそも仮説で因子数2と仮

因子負荷量・因子得点の算出

関数を定義

# 因子負荷量と因子得点を算出する関数を作成
def fa_promax(data,n,rotation):
    
    # データの標準化
    data = (data - data.mean())/data.std()
    fa = FA(n_factors=n, method='principal', rotation=rotation)
    fa.fit(data)
    
    # カラム名の作成
    columns = []
    for i in range(n):
        colnum = i + 1
        colname = '第{}因子'.format(colnum)
        columns.append(colname)
        
    fig = plt.figure(fi視聴gsize=(15,7))
    ax1 = fig.add_subplot(1,2,1)
    ax2 = fig.add_subplot(1,2,2)
    
    # 回転後の因子負荷量の表示
    loadings_df = pd.DataFrame(fa.loadings_, columns=columns) 
    loadings_df.index = data.columns

    sns.heatmap(loadings_df,annot=True, fmt='.2f',cmap='GnBu',ax=ax1)

    # 横軸に第1因子負荷量、縦軸に第2因子負荷量でプロット
    ax2.scatter(loadings_df['第1因子'],loadings_df['第2因子'])
    ax2.plot([0,0] , [-1,1] ,linestyle='--',color='grey')
    ax2.plot([-1,1] , [0,0] ,linestyle='--',color='grey')
    ax2.set_xlim(-1,1)
    ax2.set_ylim(-1,1)
    ax2.set_xlabel('第1因子負荷量')
    ax2.set_ylabel('第2因子負荷量')

    # ラベルを表示
    for i in range(len(loadings_df)):
        ax2.text(loadings_df.iloc[i][0], loadings_df.iloc[i][1], loadings_df.index[i])

    plt.show()
    
    # 因子得点を返す
    return pd.DataFrame(fa.transform(data),columns=columns)

作成した関数を実行

result = fa_promax(data = data.drop(['年齢','年齢層'],axis=1),
                   n = 2, # 因子数=2
                   rotation='promax') # promax回転

第1因子は、「コンテンツの見つけやすさ」～「視聴の質」が強く影響しており、想定通りアプリの機能性と一致する。第1因子の定義を「視聴性能への期待感」とした。
第2因子は、「コンテンツの充実」～「コンテンツの更新頻度」「価格」が強く影響しており、想定通りコンテンツの充実感と一致する。第2因子の定義を「お得感」とした。

参考までに）回転なしの場合

第1因子が全て+、第2因子に±が混在している。これでは、第1因子は「総合評価」になってしまい、第2因子は解釈しづらい。これを回転してみると、以下のとおり。

Promax回転は、斜交回転のため、回転軸が直交していない（赤い軸同士は直角でない）ことが以下からもわかると思う（ほぼ直角でわかりづらいですが、気持ち直角でないことを感じていただければ。）

ちなみにPromax回転は、一度Varimax回転（直交回転）を経由し、そこからさらに斜交回転を加える形になるので、以下のイメージの方が適切かもしれません。

年齢層別で因子分析

年齢層でデータフレームを区切り、因子分析をかけた結果のみ記載します。スクリ―プロットは記載しておりませんが、描画の結果、全体の場合と同様の判断で、因子数を2としています。

20代以下（18~29歳）

他の年齢層と比べると、第2因子に「価格」があまり貢献していないことがわかります。

30-40代（30~49歳）

50代以上（50~69歳）

他の年齢層と比べると、第2因子に「ジャンルの充実」があまり貢献していないことがわかります。

結果

①：現状のサービス満足度を評価する
50代以上は「ジャンルの充実」満足度、20代以下では「価格」満足度が低い
②：各満足度の背景にある因子を推察し、仮説と相違ないか検証する
仮説で設定した因子：「アプリの機能性」、「コンテンツの充実感」に対し、実際に定義した因子：「視聴性能への期待感」、「お得感」であるため、大きな乖離はなく、仮説どおりの因子を抽出できた
③：年齢層別で傾向が異なるか比較する
第2因子「お得感」に対して、20代以下では「価格」が、50代以上では「ジャンルの充実」があまり貢献していなかった。即ち、20代以下では「価格」に対する満足度と「お得感」は他の年齢に比べて相関が弱く、50代以上では「ジャンルの充実」に対する満足度と「お得感」は他の年齢に比べて相関が弱い。
④：①～③から今後も継続してもらうための改善示唆を出す
（少し強引ですが）若年者割引プランの策定や、高年齢層が興味あるジャンルの調査とその反映等の改善アクションによって、現ユーザーのさらなる継続利用が見込める可能性あり

補足：因子分析（Factor-analyzerなし）

ライブラリ「Factor-analyzer」を使用せずに因子分析を実行してみる（因子数2の場合）。
なお、以下①では主因子法、共通性の計算は線形回帰を使用、③ではターゲット行列の計算にκ係数=4としています。

【導出過程】
①：回転前の因子負荷量を計算
②：バリマックス回転後の因子負荷量を計算
③：プロマックス回転後の因子負荷量を計算

①：回転前の因子負荷量を計算

# 因子数=2とした場合の算出関数
def fa_nonrotation(data):
    
    # データの標準化
    data = (data - data.mean())/data.std()
    
    # 相関行列の作成
    corr = data.corr()
    
    # 共通性（1-d^2）の算出・・・線形回帰を採用
    cols = data.columns
    for col in cols:
        model = sm.OLS(data[col],data.drop(col,axis=1)).fit()
        val = model.rsquared
        
        # 相関行列の値を置換（自分対自分）
        corr[col][col] = val
    
    # 固有値、固有ベクトルを求める
    eig_val, eig_vec = np.linalg.eig(corr)
    
    # 固有値に対応する固有ベクトルを辞書に格納
    eig_dict = {}
    for i in range(len(eig_val)):
        key = eig_val[i]
        vec_i = eig_vec[:,i]
        eig_dict[key] = vec_i
    
    # 固有値を降順で並べ替え、因子数分を頭から取得
    eig_val = list(eig_dict.keys())
    eig_val.sort(reverse=True)
    eig_val = np.array(eig_val)
    eig_val = eig_val[0:2]
    
    # 上記で取得した固有値に対応する固有ベクトルを取得
    vec_vals = []
    for i in eig_val:
        vec_vals.append(list(eig_dict[i]))
    eig_vec = np.array(vec_vals)
    
    # 因子負荷量から相関係数を逆算する式で戻す（行列の内積ではなく、同一項の単純な積）
    A1 = eig_vec.T*np.sqrt(eig_val)
    A2 = (eig_vec.T*np.sqrt(eig_val)).T
    
    # A1とA2の内積を計算
    dot = np.dot(A1,A2)
    
    # 0回目で超過はないと思うが判定を実施
    for i in range(len(dot)):
        if dot[i][i]>1:
            print('0回目で任意の共通性が1を超えました')
            return
            
        else:
            pass
            
    # 共通性の値とA1とA2の内積の対角値を置き換え
    for i in range(len(corr)):
        corr.iloc[i][i] = dot[i][i]

    # 因子負荷量を辞書型に保存していく
    fa_result = {}
    counter = 0
    fa_result[counter] = A1

    #----------上記までの処理で対角値>1はない----------
    # ここからループ処理で、対角値（任意の共通性）>1が出るまで続け、その1つ前の因子負荷量を採用する
    loop_count = 50
    for c in range(loop_count):

        # 固有値、固有ベクトルを求める
        eig_val, eig_vec = np.linalg.eig(corr)

        # 固有値に対応する固有ベクトルを辞書に格納
        eig_dict = {}
        for i in range(len(eig_val)):
            key = eig_val[i]
            vec_i = eig_vec[:,i]
            eig_dict[key] = vec_i

        # 固有値を降順で並べ替え、因子数分を頭から取得
        eig_val = list(eig_dict.keys())
        eig_val.sort(reverse=True)
        eig_val = np.array(eig_val)
        eig_val = eig_val[0:2]

        # 上記で取得した固有値に対応する固有ベクトルを取得
        vec_vals = []
        for i in eig_val:
            vec_vals.append(list(eig_dict[i]))
        eig_vec = np.array(vec_vals)

        # 因子負荷量から相関係数を逆算する式で戻す（行列の内積ではなく、同一項の単純な積）
        A1 = eig_vec.T*np.sqrt(eig_val)
        A2 = (eig_vec.T*np.sqrt(eig_val)).T

        # A1とA2の内積を計算
        dot = np.dot(A1,A2)

        # 共通性>1の判定を実施
        comon_counter = 0
        for i in range(len(dot)):
            if dot[i][i]>1:
                print('{}回目で任意の共通性が1を超えました'.format(counter+1))
                comon_counter += 1
                
                # 辞書に格納されているキーの中で最大のものを抽出
                key_list = list(fa_result.keys())
                key_list.sort(reverse=True)
                key = key_list[0]
                
                # 辞書中の値を取得
                val = fa_result[key]
                result = pd.DataFrame(val,columns=['第1因子', '第2因子'],index=data.columns)
                
                
                fig = plt.figure(figsize=(15,7))
                ax1 = fig.add_subplot(1,2,1)
                ax2 = fig.add_subplot(1,2,2)
                sns.heatmap(result,annot=True, fmt='.2f',cmap='GnBu',ax=ax1)
                
                # 横軸に第1因子負荷量、縦軸に第2因子負荷量でプロット
                ax2.scatter(result['第1因子'],result['第2因子'])
                ax2.plot([0,0] , [-1,1] ,linestyle='--',color='grey')
                ax2.plot([-1,1] , [0,0] ,linestyle='--',color='grey')
                ax2.set_xlim(-1,1)
                ax2.set_ylim(-1,1)
                ax2.set_xlabel('第1因子負荷量')
                ax2.set_ylabel('第2因子負荷量')
                
                # ラベルを表示
                for i in range(len(result)):
                    ax2.text(result.iloc[i][0], result.iloc[i][1], result.index[i] )
                
                plt.show()
                
                # 関数の返り値は因子負荷量
                return result

            else:
                pass
        
        # 共通性の値とA1とA2の内積の対角値を置き換え
        for i in range(len(corr)):
            corr.iloc[i][i] = dot[i][i]

        # 因子負荷量を辞書型に保存していく
        counter += 1
        fa_result[counter] = A1
        
    # どれだけ繰り返しても共通性が1を超えなかった場合（最終処理を採用する）
    if comon_counter == 0:
        print('{}回の試行では任意の共通性は1を超えませんでした'.format(loop_count))
        
        # 辞書に格納されているキーの中で最大のものを抽出
        key_list = list(fa_result.keys())
        key_list.sort(reverse=True)
        key = key_list[0]

        # 辞書中の値を取得
        val = fa_result[key]
        result = pd.DataFrame(val,columns=['第1因子', '第2因子'],index=data.columns)
        
        fig = plt.figure(figsize=(15,7))
        ax1 = fig.add_subplot(1,2,1)
        ax2 = fig.add_subplot(1,2,2)
        sns.heatmap(result,annot=True, fmt='.2f',cmap='GnBu',ax=ax1)

        # 横軸に第1因子負荷量、縦軸に第2因子負荷量でプロット
        ax2.scatter(result['第1因子'],result['第2因子'])
        ax2.plot([0,0] , [-1,1] ,linestyle='--',color='grey')
        ax2.plot([-1,1] , [0,0] ,linestyle='--',color='grey')
        ax2.set_xlim(-1,1)
        ax2.set_ylim(-1,1)
        ax2.set_xlabel('第1因子負荷量')
        ax2.set_ylabel('第2因子負荷量')

        # ラベルを表示
        for i in range(len(result)):
            plt.text(result.iloc[i][0], result.iloc[i][1], result.index[i] )

        plt.show()
        return result
        
    else:
        pass

# 関数を実行
result_df = fa_nonrotation(data=data.drop(['年齢','年齢層'],axis=1))

②：バリマックス回転後の因子負荷量を計算

②-1：回転角度の決定

# 回転前の因子負荷行列をA1に格納
A1 = result_df.copy()
A1 = np.array(A1)

# Varimaxによる最大化する値の格納リストを作成
varimax_list = []

# 回転するΘは-90~0°で設定
for theta in range(-180,181,1):
    
    sin = math.sin(math.radians(theta))
    cos = math.cos(math.radians(theta))
    angle = np.array([[cos,sin],
                      [-sin,cos]])
    
    # 回転前の因子負荷行列と回転行列の内積を計算
    B1= np.dot(A1,angle)
    
    # 最大化する値の計算のため、行列値を二乗
    B1 = B1**2
    B1 = pd.DataFrame(B1)
    B1['0+1'] = B1[0] + B1[1]
    B1['0/0+1'] = B1[0]/B1['0+1']
    B1['1/0+1'] = B1[1]/B1['0+1']

    # 合計値計算用に空のリストを作成
    nums = []
    for i in range(len(B1)): 
        num1 = ( B1[0][i] / (B1[0][i] + B1[1][i]) ) - sum(B1['0/0+1'])/len(B1)
        num1 = num1**2
        num2 = ( B1[1][i] / (B1[0][i] + B1[1][i]) ) - sum(B1['1/0+1'])/len(B1)
        num2 = num2**2
        num = num1 + num2
        nums.append(num)
    sum_num = sum(nums)
    add_num = [theta,sum_num]
    varimax_list.append(add_num)
    
varimax_results = pd.DataFrame(varimax_list,columns=['theta','varimax'])
display(varimax_results[varimax_results['varimax']==varimax_results['varimax'].max()])

回転角度-107°が見やすいので採用

②-2：回転後の因子負荷量

# varimax回転後の因子負荷行列を算出
theta = -107
sin = math.sin(math.radians(theta))
cos = math.cos(math.radians(theta))
angle = np.array([[cos,sin],
                  [-sin,cos]])

B1= np.dot(A1,angle)
result = pd.DataFrame(B1,columns=['第1因子', '第2因子'],index=data.drop(['年齢','年齢層'],axis=1).columns)

     
fig = plt.figure(figsize=(15,7))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)

sns.heatmap(result,annot=True, fmt='.2f',cmap='GnBu',ax=ax1)

# 横軸に第1因子負荷量、縦軸に第2因子負荷量でプロット
ax2.scatter(result['第1因子'],result['第2因子'])
ax2.plot([0,0] , [-1,1] ,linestyle='--',color='grey')
ax2.plot([-1,1] , [0,0] ,linestyle='--',color='grey')
ax2.set_xlim(-1,1)
ax2.set_ylim(-1,1)
ax2.set_xlabel('第1因子負荷量')
ax2.set_ylabel('第2因子負荷量')

# ラベルを表示
for i in range(len(result)):
    ax2.text(result.iloc[i][0], result.iloc[i][1], result.index[i])

plt.show()

③：プロマックス回転後の因子負荷量を計算

③-1：ターゲット行列の計算

ターゲット行列：
本来の因子負荷量はこれくらいだろう、とする行列。プロマックス回転では、このターゲット行列にできるだけ近くなるように値を変換している。

# ターゲット行列の計算
result_tmp = result.copy()
result_tmp = np.array(result_tmp)
result_tmp = pd.DataFrame(result_tmp)
result_tmp['√0+1'] = np.sqrt(result_tmp[0]**2 + result_tmp[1]**2)
target = []
K=4 # カッパ係数
for i in range(len(result_tmp)):
    num0 = ( result_tmp['√0+1'][i]/result_tmp[0][i] ) * ((np.abs( result_tmp[0][i]/result_tmp['√0+1'][i] ))**(K+1))
    num1 = ( result_tmp['√0+1'][i]/result_tmp[1][i] ) * ((np.abs( result_tmp[1][i]/result_tmp['√0+1'][i] ))**(K+1))
    nums = [num0,num1]
    target.append(nums)
    
print('ターゲット行列')
display(pd.DataFrame(target,columns=['第1因子','第2因子'],index=result.index).round(5))

③-2：回転後の因子負荷量

# バリマックスの因子負荷量行列
aB = np.array(result)

# ターゲット行列
aTg = target

# 行列Q,d,D
aQ = np.dot(np.linalg.inv(np.dot(aB.T,aB)),np.dot(aB.T,aTg))
ad = np.linalg.inv(np.dot(aQ.T,aQ))
aD = np.diag(np.diag(ad, k=0))
aD = np.sqrt(aD)

# プロマックスの因子負荷量行列
aP = np.dot(np.dot(aB,aQ),aD)
print('プロマックスの因子負荷量行列')
aP = pd.DataFrame(aP,columns=['第1因子','第2因子'],index=result.index)

fig = plt.figure(figsize=(15,7))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)

sns.heatmap(aP.round(5),annot=True, fmt='.2f',cmap='GnBu',ax=ax1)

# 横軸に第1因子負荷量、縦軸に第2因子負荷量でプロット
ax2.scatter(aP['第1因子'],aP['第2因子'])
ax2.plot([0,0] , [-1,1] ,linestyle='--',color='grey')
ax2.plot([-1,1] , [0,0] ,linestyle='--',color='grey')
ax2.set_xlim(-1,1)
ax2.set_ylim(-1,1)
ax2.set_xlabel('第1因子負荷量')
ax2.set_ylabel('第2因子負荷量')

# ラベルを表示
for i in range(len(aP)):
    ax2.text(aP.iloc[i][0], aP.iloc[i][1], aP.index[i])


plt.show()

Factor-analyzer使用時と比較して、因子が逆になってはいることと、因子負荷量が全く同じというわけにはいきませんでしたが、近しい値になっていると思います。