「入門はじめての多変量解析」をPythonで写経 Vol.16 ～ ４章「はじめての因子分析」①Hello, 因子分析！

４章「はじめての因子分析」

書籍の著者　石村貞夫 先生、石村光資郎 先生

書籍「入門はじめての多変量解析」４章「はじめての因子分析」の Python写経活動記録 です。 

多変量解析の入門を Python と一緒に学ぶ写経シリーズです。

この記事は、仮想データで ざっくり因子分析 に取り組みます。
書籍の因子分析に突入する前のウォーミング・アップの位置づけです。
因子分析のうち「探索的因子分析」を取り扱います。

ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

４章 はじめての因子分析

この記事で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd

# 因子分析
from factor_analyzer import FactorAnalyzer
from factor_analyzer.factor_analyzer import (
    calculate_kmo, calculate_bartlett_sphericity)

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

# ワーニングの非表示 ※factor_analyzer実行時に発生
import warnings
warnings.filterwarnings('ignore')

因子分析の一連の流れを体感する

この記事は、架空のアンケート調査結果で因子分析を行い、データの背後にある「共通因子」（潜在的な変数）を探ります。

ChatGPTによると：

因子分析は、観測された複数の変数どうしの「共通して動いている部分（共分散／相関）」に着目し、そこから

1. 背後にある潜在的な要素（共通因子）が何本あるか

2. 各変数がそれらの因子のどれに、どのくらい影響を受けているか

を明らかにする手法です。

データの分散共分散・相関係数に着目して、データの変数が影響を受けている「共通因子（因子と略すことも）」を見つけに行きましょう！

１．データの概要
ChatGPTと相談しながら作成した仮想のアンケートデータです。

### 1. データサンプル生成

# 設定
n_samples = 50
rng = np.random.default_rng(seed=2) # 2, 3, 7

# 潜在因子を3つ生成
f1 = rng.normal(size=n_samples)
f2 = rng.normal(size=n_samples)
f3 = rng.normal(size=n_samples)

# 各観測変数への負荷量行列
# item1-3 -> f1, item4-5 -> f2, item6-8 -> f3
load_matrix = np.array([
    [0.8, 0, 0],  # Q1 味
    [0.7, 0, 0],  # Q2 香り
    [0.9, 0, 0],  # Q3 ブランド
    [0, 0.7, 0],  # Q4 量
    [0, 0.8, 0],  # Q5 値段
    [0, 0, 0.7],  # Q6 ストーリー性
    [0, 0, 0.8],  # Q7 飲む時間
    [0, 0 ,0.6]   # Q8 飲む場所
])

# 観測データX = 潜在因子@負荷量 + ノイズ
factors = np.column_stack([f1, f2, f3])  # (n_samples,3)
noise = rng.normal(scale=0.3, size=(n_samples, 8))
X = factors @ load_matrix.T + noise

# DataFrame化・5段階評価に丸め込み
index = [f'No.{i+1:02}' for i in range(n_samples)]
columns = ['味', '香り', 'ブランド', '量', '値段', 'ストーリー', '時間', '場所']
df = pd.DataFrame(X, columns=columns, index=index
                  ).apply(lambda x: np.clip(np.round(x + 3), 1, 5).astype(int))

# 結果の表示
print('df.shape:', df.shape)
df.head(10)

【実行結果】
「コーヒーのこだわり」に関するアンケートを50人に行った結果、と見立てましょう。
行が回答者番号です。標本サイズは $${n=50}$$です。
列が質問項目です。項目数は $${p=8}$$です。

■ アンケートの概要
「コーヒーの◯◯について、どの程度こだわりを持っていますか」の設問に対して、以下の５の選択肢から１つを選択して回答するものです。
　１：こだわらない
　２：ややこだわらない
　３：どちらとも言えない
　４：ややこだわる
　５：こだわる

「コーヒーの◯◯について」は具体的には以下の内容です。

回答は例えば「時間が４」の場合、コーヒーを飲む時間（とき）に「ややこだわる」です。

２．データの相関行列の確認
ChatGPTによると

データの相関行列に着目しましょう。

### 2. 相関行列の確認
corr = df.corr()
print('データの相関行列:')
corr.round(3)

【実行結果】

相関の強さでなんとなく変数のグループができそうです。

• 「味・香り・ブランド」：互いに中程度以上の相関

• 「量・値段」：互いに中程度の相関

• 「ストーリー・時間・場所」：互いに中程度以上の相関

３．データの固有値とスクリープロット
因子分析の準備に入ります。
分析前に「因子の数」にあたりをつけましょう。

ここでは
　・データの「相関行列の固有値」の大きさ
　・スクリープロット
で数を決めてみます。

### 3. 固有値の計算とスクリープロット

## データの相関行列の固有値を算出
eigvals, _ = np.linalg.eig(corr)
eigvals_sorted = np.sort(eigvals)[::-1]  # 降順ソート
print('固有値:\n', np.round(eigvals_sorted, 3))

## スクリープロットを描画
# 描画領域の設定
plt.figure(figsize=(6, 4))
# 固有値の折れ線グラフの描画
plt.plot(range(1, len(eigvals_sorted)+1), eigvals_sorted, marker='o')
# 固有値=1の水平線の描画
plt.axhline(1, color='tab:red', linestyle='--')
# 修飾
plt.title('スクリープロット')
plt.xlabel('因子番号', fontsize=12)
plt.ylabel('固有値', fontsize=12)
plt.show()

【実行結果】
因子１～３と４以降で明らかに何かが違っている感じがしますね！

因子の数を「３」にして分析を進めます！

４．因子分析の実行、因子負荷量の算出
早速、因子分析を実行しましょう。
factor_analyzer ライブラリの FactorAnalyzer() を利用します。
あわせて因子負荷量を算出します。
詳しいことは後ほど。

### 4. 因子分析の実行：因子数３、主因法（テキストの最初の推定法）、バリマックス回転

# 因子分析の実行
fa = FactorAnalyzer(n_factors=3, method='principal', rotation='varimax')
fa.fit(df)

# 因子負荷量の算出
loadings = pd.DataFrame(
    fa.loadings_, index=columns, columns=['因子1','因子2','因子3'])
print('因子負荷量:')
loadings.round(3)

【実行結果】
因子負荷量は次の値になりました。

少々込み入った内容になりますが、因子分析の中身を確認します。

■ 因子分析の概要
FactorAnalyzer の引数に次の３つを設定しています。

■ 数式で見る因子分析
因子分析は重回帰モデルのような計算式で表すことができます。
ざっくり、「変数を共通因子と独自因子で説明するモデル」です。

数式で書くと…
データが標準化されているものとして、$${n}$$ 個のサンプルの $${i}$$ 番目、$${p}$$ 個の変数の $${j}$$ 番目の変数の値 $${x_{ij}}$$について、次のように示されます。

$${f_{ik}}$$ は $${m}$$個の共通因子 $${f_{k}}$$ のサンプル $${i}$$ に対する因子得点、$${a_{jk}}$$ は共通因子 $${f_{k}}$$ の変数 $${x_j}$$ に対する因子負荷量、$${\varepsilon_{ij}}$$ はサンプル $${i}$$ の変数 $${x_j}$$ に対する独自因子 です。

テキストでは、p.128 の「その３」で「変数の数が $${p}$$ 個」、「共通因子の数が $${m}$$ 個」の場合のモデル式を次のように紹介しています。
$${x_j, f_k, \varepsilon_j}$$ は $${n}$$ 次元のベクトルと思われます。

■ 因子負荷量 $${a_{jk}}$$ の推定
上の数式を参考にすると…
右辺の「因子負荷量」「因子得点」「独自因子」の全ては、決まった値が与えられていない変数です。
そのため、因子負荷量を推定するには他の変数の推定も必要になり、推定計算が複雑になります。

そのような推定方法には複数の方法があるようです。

• テキストは「主因法」と「最尤法」を紹介しています。

• FactorAnalyzer() は「MINRES法」（最小残差法：デフォルト）、主因法、最尤法に対応しています。

• statsmodels には因子分析 multivariate.factor.Factor() があり、主因法（デフォルト）と最尤法に対応しています。

• scikit-learn には因子分析 FactorAnalysis() があり、最尤法に対応しています。

■ 回転
因子分析の特徴の１つが「回転」です。
因子負荷量を文字通り回転するものです。
回転により「共通因子に強く影響する変数」が把握しやすくなります。
回転後は「各変数ができるだけ１つの因子に強く関連し、それ以外の因子はほぼゼロに近づく」からです。
結果として、共通因子の解釈がしやすくなります。

■ 共通因子の解釈
先ほど Python で表示した因子負荷量は回転済みの状態です。
解釈がしやすくなったのか、次のコードで確かめてみましょう。

# 因子負荷量の可視化
sns.heatmap(loadings, annot=True, fmt='.3f', cmap='Greens',
            annot_kws=dict(fontsize=12))
plt.title('因子負荷量');

【実行結果】
列に因子番号１～３が並び、行に変数が並んでいます。

共通因子ごとに因子負荷量が大きくなっている変数が見つかります！
共通因子１は「味・香り・ブランド」に強い影響を与えています。
共通因子２は「ストーリー・時間・場所」に、共通因子３は「量・値段」に強い影響を与えています。

変数の意味合いを想像しながら、共通因子の解釈を試みてみましょう。

【解釈（例）】

• 因子1 (F1)：味・香り・ブランド に強くロード
  → 共通因子１はコーヒーの「クオリティ」

• 因子2 (F2)：ストーリー・時間・場所 に強くロード
  → 共通因子２はコーヒーを楽しむ「シーン」

• 因子3 (F3)：量・値段 に強くロード
  →共通因子３はコーヒーに求める「コスパ」

みなさんはどんな解釈をしますか？

５．因子得点の算出
FactorAnalyzer で回答者ごとの因子得点 $${f_{ki}}$$ を算出しましょう。

### 5. 因子得点の算出 ※FactorAnalyzerの推定法は「バートレットの方法」とのこと
scores = pd.DataFrame(fa.transform(df), index=df.index,
                      columns=['因子得点1','因子得点2','因子得点3'])
print('因子得点:')
scores.head(10).round(3)

【実行結果】
先頭 10 名の因子得点を表示しています。

• 因子得点１が高い人ほど「クオリティの因子」を強く持つ

• 因子得点２が高い人ほど「シーンの因子」を強く持つ

• 因子得点３が高い人ほど「コスパの因子」を強く持つ

という解釈ができます。

６．バイプロットの描画
回答者ごとの因子得点と変数ごとの因子負荷量を重ねてプロットします。

こちらは共通因子１と共通因子２の描画です。

### 6.バイプロットの描画 F1, F2

# 描画領域の設定
fig, ax = plt.subplots(figsize=(7, 7))
twin = ax.twinx().twiny()

# 因子得点のテキストの描画
for x, y, s in zip(scores['因子得点1'], scores['因子得点2'], scores.index):
    ax.text(x=x, y=y, s=s, ha='center', va='center', color='tab:blue')
# x=0の垂直線、y=0の水平線の描画
ax.axvline(0, color='black', lw=0.5)
ax.axhline(0, color='black', lw=0.5)
# 修飾
ax.set(xlim=(-3, 3), ylim=(-3, 3))
ax.set_xlabel('第１因子: クオリティ', fontsize=12)
ax.set_ylabel('第２因子: シーン', fontsize=12)

# 因子負荷量の矢印と質問名の描画
for x, y, text in zip(loadings['因子1'], loadings['因子2'], loadings.index):
    # annotateで矢印をプロット
    twin.annotate(
        text='',          # テキスト：index
        xy=[x, y],        # テキストの位置:[第1因子負荷量, 第2因子負荷量]
        xytext=[0, 0],    # 終点の位置：[0, 0]
        color='tab:red',  # 文字の色
        arrowprops=dict(  # 矢印の設定
            arrowstyle='->',       # 矢印の形状と向き（逆向きに設定しています）
            facecolor='tab:red',   # 矢印の表面の色
            edgecolor='tab:red',   # 矢印の輪郭の色
        )
    )
    # testでテキストを描画
    twin.text(x=x, y=y, s=text, color='tab:red', fontsize=12)

# 修飾
twin.set(xlim=(-1.2, 1.2), ylim=(-1.2, 1.2))
plt.show()

【実行結果】
「クオリティ」に関しては、チャートの右側が重要視する人、左側が重要視しない人です。
コーヒーを飲む「シーン」に関しては、チャートの上側が重要視する人、下側が重要視しない人です。

こちらは共通因子１と共通因子３の描画です。

### 6.バイプロットの描画 F1, F3

# 描画領域の設定
fig, ax = plt.subplots(figsize=(7, 7))
twin = ax.twinx().twiny()

# 因子得点のテキストの描画
for x, y, s in zip(scores['因子得点1'], scores['因子得点3'], scores.index):
    ax.text(x=x, y=y, s=s, ha='center', va='center', color='tab:blue')
# x=0の垂直線、y=0の水平線の描画
ax.axvline(0, color='black', lw=0.5)
ax.axhline(0, color='black', lw=0.5)
# 修飾
ax.set(xlim=(-3, 3), ylim=(-3, 3))
ax.set_xlabel('第１因子: クオリティ', fontsize=12)
ax.set_ylabel('第３因子: コスパ', fontsize=12)

# 因子負荷量の矢印と質問名の描画
for x, y, text in zip(loadings['因子1'], loadings['因子3'], loadings.index):
    # annotateで矢印をプロット
    twin.annotate(
        text='',          # テキスト：index
        xy=[x, y],        # テキストの位置:[第1因子負荷量, 第2因子負荷量]
        xytext=[0, 0],    # 終点の位置：[0, 0]
        color='tab:red',  # 文字の色
        arrowprops=dict(  # 矢印の設定
            arrowstyle='->',       # 矢印の形状と向き（逆向きに設定しています）
            facecolor='tab:red',   # 矢印の表面の色
            edgecolor='tab:red',   # 矢印の輪郭の色
        )
    )
    # testでテキストを描画
    twin.text(x=x, y=y, s=text, color='tab:red', fontsize=12)

# 修飾
twin.set(xlim=(-1.2, 1.2), ylim=(-1.2, 1.2))
plt.show()

【実行結果】
「コスパ」に関しては、チャートの上側が重要視する人、下側が重要視しない人です。

仮想データは共通因子を見つけやすいように「意図的に」構築しています。
では日常のデータで因子分析を実践するときに、どんなデータが因子分析に適していて、好ましい因子分析の結果はどのようなものなのでしょう？

因子分析を評価する

因子分析に関するさまざまな評価をChatGPTに訊きました。
評価ポイントの概要、仮想データを使った評価例を確認しましょう。

【注意点】
ChatGPTが回答した 評価手法の妥当性・網羅性や指標の目安の適否は未確認 ですので、ご留意ください。

最初にデータそのものの評価です。

１．KMO（Kaiser–Meyer–Olkin）検定

• 何を測る？

  • 項目間の相関に対する「部分相関」の大きさを比べて、因子構造に適したデータかどうかを示す指標。

• 値の目安

  • 0.90～1.00：  極めて良好

  • 0.80～0.89：  非常に良い

  • 0.70～0.79：  良い

  • 0.60～0.69：  可（最低限の許容水準）

  • 0.50～0.59：  不十分（再検討）

  • 0.00～0.49：  不適

### 1. KMO テスト
kmo_per_item, kmo_overall = calculate_kmo(df)
print(f'Overall KMO: {kmo_overall:.3f}')
print(f'KMO per item min: {kmo_per_item.min():.3f}')
pd.Series(kmo_per_item, index=columns, name='KMO').to_frame().round(3)

【実行結果】
全体の KMO は「可」です。よしとしましょう！

２．Bartlett の球面性検定

• 何を測る？

  • 「相関行列が単位行列（互いに無相関）である」という帰無仮説を検定。

• 値の目安

  • $${p < 0.05}$$： 因子分析に必要な相関がある → 適

  • $${p \geq 0.05}$$： 相関が弱く因子構造が乏しい → 不適

### 2. Bartlett の球面性検定
chi_square_value, p_value = calculate_bartlett_sphericity(df)
print(f"Bartlett's Test: χ² = {chi_square_value:.1f}, p = {p_value:.3g}\n")

【実行結果】
$${p}$$ 値は $${<0.05}$$ です。
有意水準 $${5\%}$$ で有意であり、因子分析に必要な相関があると言えそうです。

ここからは共通因子の数の判断材料です。
既出の項目も改めて書いています。

３．固有値（データの相関行列に関する）

• 何を測る？

  • 抽出した各因子が「標本データの総分散（情報量）をどれだけ持っているか」。

• 値の目安

  • $${λ > 1}$$： その因子は「１変数分以上の情報量」をもつ
    → 採用の目安（Kaiser基準）

  • $${λ < 1}$$： 「１変数分以下」の情報しかない
    → 棄却

### 3. データの相関行列の固有値
eigvals, _ = np.linalg.eig(corr)
eigvals_sorted = np.sort(eigvals)[::-1]
print('データの相関行列の固有値:')  # スクリープロットで使用する固有値
print(eigvals_sorted.round(3))

【実行結果】
$${λ > 1}$$ の固有値の個数は $${3}$$ です。
因子の数は $${3}$$ を推します！

４．スクリープロット（Scree Plot）

• 何を測る？

  • 固有値を降順に折れ線グラフ化し、「肘（エルボー）」で因子数を判断。

• 目安

  • 折れ線が急に水平に近づく「肘の位置」の直前までを採用

### 4. スクリープロット

# 描画領域の設定
plt.figure(figsize=(6, 4))
# 固有値の折れ線グラフの描画
plt.plot(range(1, len(eigvals_sorted)+1), eigvals_sorted, marker='o')
# 固有値=1の水平線の描画
plt.axhline(1, color='tab:red', linestyle='--')
# 修飾
plt.title('スクリープロット')
plt.xlabel('因子番号', fontsize=12)
plt.ylabel('固有値', fontsize=12)
plt.show()

【実行結果】
肘の部分は $${4}$$ です（綺麗に出過ぎ！？）。
因子の数は $${3}$$ を推します！

５．平行分析（Parallel Analysis）

• 何を測る？

  • 実データの固有値と、同じサイズのランダムデータを何度も生成して得られる固有値を比較し、「偶然の相関で得られる固有値」と「実データの固有値」を区別 します。

  • 具体的にはランダムデータの固有値の分布（たとえば 95th パーセンタイル）を基準にし、実データの固有値がその閾値を超える主成分だけを「真の因子」とみなします。

• 値の目安

  • ランダム閾値：ランダムデータの固有値分布の 95th パーセンタイルをよく使う

  • 因子採用の条件：
    「$${\text{実データの固有値}_i > \text{ランダムデータの95\%点の固有値}_i}$$」を満たす最初の $${i}$$ までを因子数として採用

  • 繰り返し数（n_iter）：最低 100～200、できれば 500～1,000 回以上で安定化

これにより、固有値 ＞ 1 ルールやスクリープロットの主観を補い、「偶然のノイズ以上に有意な因子」を客観的に選定できます。

平行分析を実行できる適当なライブラリが見当たらなかったので、ChatGPTと一緒に関数化に取り組みました。

【方法１】主成分を対象にした平行分析
固有値には「データの相関行列の固有値」を用いています。

### 5. 平行分析

## 主成分を対象にした平行分析関数の定義
# データの相関係数の固有値を使う。因子分析後のFAの固有値ではない

def parallel_analysis(df, n_iter=500, quantile=0.95, seed=None):
    
    ## 設定と準備
    # 乱数生成器
    rng = np.random.default_rng(seed)
    # データの標本サイズ、説明変数の数
    n, p = df.shape
    
    ## データの相関行列の固有値の算出
    real_eig, _ = np.linalg.eig(df.corr())
    real_eig = np.sort(real_eig)[::-1]  # 降順ソート

    ## ランダム行列の相関行列の固有値の算出
    # 結果を格納する配列の準備
    rand_eigs = np.zeros((n_iter, p))
    # ランダム行列を生成してi_iter個の固有値を算出
    for i in range(n_iter):
        # ランダム行列データ：標準正規分布乱数
        rand = rng.standard_normal(size=(n, p))
        # 相関行列の固有値の算出
        eigs, _ = np.linalg.eig(np.corrcoef(rand, rowvar=False))
        rand_eigs[i, :] = np.sort(eigs)[::-1] # 降順ソート

    ## 結果のまとめ
    # ランダム行列の固有値のパーセンタイルをしきい値とする
    thresh = np.quantile(rand_eigs, quantile, axis=0)
    # データの固有値とランダム行列の固有値のパーセンタイルをデータフレーム化
    pa = pd.DataFrame({
        'Real PC EigenValue': real_eig, f'Random PC Q{quantile}': thresh
        }, index=[f'因子{i+1}' for i in range(p)])
    # 平行分析の推奨因子数の算出 ※ランダム行列の95パーセンタイルよりも大きい個数
    suggested_factors = sum(real_eig > thresh)
    
    ## 戻り値：固有値のデータフレーム、推奨因子数
    return pa, suggested_factors

【実行結果】なし

平行分析を実行します。

## 主成分を対象にした平行分析による因子数の算出
result = parallel_analysis(df, n_iter=500, seed=0)
print('平行分析による推奨因子数 =', result[1])

【実行結果】
推奨は $${3}$$ です！

こちらもスクリープロットで可視化できます。

## スクリープロットの描画

# 描画領域の設定
plt.figure(figsize=(6, 4))
# データの固有値とランダム行列の95パーセンタイルの折れ線グラフの描画
sns.lineplot(result[0], marker='o', ms=8)
# 固有値=1の水平線の描画
plt.axhline(1, color='gray', lw=1, ls='--')
# 修飾
plt.title('スクリープロット\nデータのPC固有値, ランダム行列のPC固有値の95%点')
plt.ylabel('固有値', fontsize=12)
plt.show()

【実行結果】
「ランダムを示すオレンジの点線」より上にある「実データの青い点」は３つです！

【方法２】因子を対象にした平行分析
固有値には「因子の固有値」を用いています。

## 因子を対象にした平行分析を行う関数の定義
# 因子分析後のFAの固有値を使う

def parallel_analysis_fa(df, n_iter=500, quantile=0.95, method='minres',
                         seed=None):
    """
    FA の固有値で Horn の平行分析を行う。
    
    Parameters
    ----------
    df : pandas.DataFrame
        因子分析対象のデータ（数値列のみ、標準化済みでも可）
    n_iter : int
        ランダムデータ試行回数（デフォルト500）
    quantile : float
        ランダム固有値のしきい値（デフォルト0.95）
    random_state : int or None
        乱数シード
    
    Returns
    -------
    pandas.DataFrame
        実データの FA 固有値とランダムデータのしきい値固有値
    int
        推奨因子数
    """
    
    ## 0) 設定と準備
    # 乱数生成器
    rng = np.random.default_rng(seed)
    # データの標本サイズ、説明変数の数
    n, p = df.shape

    ## 1) 実データの FA 固有値の算出
    # FactorAnalyzerのインスタンス生成
    fa_real = FactorAnalyzer(n_factors=1, rotation=None, method=method)
    # データの因子分析の実行
    fa_real.fit(df)
    # FA固有値の取得
    _, real_fa_eig = fa_real.get_eigenvalues()
    real_fa_eig = np.sort(real_fa_eig)[::-1]  # 降順ソート

    ## 2) ランダム行列で FA を繰り返し、固有値を収集
    # 結果を格納する配列の準備
    rand_eigs = np.zeros((n_iter, p))
    # ランダム行列を生成してi_iter個のFA固有値を算出
    for i in range(n_iter):
        # ランダム行列データ：標準正規分布乱数
        rand = rng.standard_normal(size=(n, p))
        # FactorAnalyzerのインスタンス生成
        fa_rand = FactorAnalyzer(n_factors=1, rotation=None, method=method)
        # ランダム行列の因子分析の実行
        fa_rand.fit(rand)
        # ランダム行列のFA固有値の取得
        _, eig_rand = fa_rand.get_eigenvalues()
        rand_eigs[i, :] = np.sort(eig_rand)[::-1]  # 降順ソート

    ## 3) 結果のまとめ
    # 3-1) ランダム行列のFA固有値を基にしてしきい値を計算（例：95パーセンタイル）
    thresh = np.quantile(rand_eigs, quantile, axis=0)
    # 3-2) 結果を DataFrame で返す
    pa = pd.DataFrame({
        'Real FA EigenValue': real_fa_eig,
        f'Random FA Q{quantile}': thresh
    }, index=[f'因子{i+1}' for i in range(p)])
    # 3-3) 平行分析の推奨因子数の算出
    suggested_factors = sum(real_fa_eig > thresh)

    ## 戻り値：FA固有値のデータフレーム、推奨因子数
    return pa, suggested_factors

【実行結果】なし

平行分析を実行します。

## 因子を対象にした平行分析による因子数の算出
result_fa = parallel_analysis_fa(df, n_iter=500, quantile=0.95, seed=0)
print('平行分析による推奨因子数 =', result_fa[1])

【実行結果】
こちらも推奨は $${3}$$ です！

スクリープロットで可視化します。

## スクリープロットの描画

# 描画領域の設定
plt.figure(figsize=(6, 4))
# データの固有値とランダム行列の95パーセンタイルの折れ線グラフの描画
sns.lineplot(result_fa[0], marker='o', ms=8)
# 修飾
plt.title('スクリープロット\nデータのFA固有値, ランダム行列のFA固有値の95%点')
plt.ylabel('固有値', fontsize=12)
plt.show()

【実行結果】
こちらも実データの最初の３点がランダムよりも上に位置しています。

【謝辞】
すうがくぶんか社の内場さま（Xアカウント @utaka233）より、平行分析の具体的な方法に関するアドバイスをいただきました。

特に、以下の点を教えていただきました。

• 主成分を対象にした平行分析と因子を対象にした平行分析が別物であることと、それらの相違点

• 因子分析を対象にした平行分析のPythonコード

いつもありがとうございます！

ここからは推定した因子負荷量の評価です。

６．共通性（Communalities）と独自性（Uniquenesses）

• 何を測る？

  • 共通性：変数ごとに「抽出した因子が説明できる分散の割合」

  • 独自性：1－共通性＝「因子が説明できなかった分散の割合」

• 値の目安

  • 共通性 ≥ 0.50： 因子が半分以上の分散を説明 → 良好

  • 共通性 < 0.40： 説明が不十分 → 要再検討

  • 独自性 < 0.50： 因子モデルとの適合が良い

### 6. 共通性・独自性　※共通性＋独自性＝１

# 共通性の取得
communalities = fa.get_communalities()
# 独自性の取得
uniquenesses = fa.get_uniquenesses()
# 結果の表示
print('共通性と独自性:')
com_uni = pd.DataFrame(
    {'共通性': communalities, '独自性': uniquenesses}, index=columns
).round(3)

com_uni

【実行結果】
共通性は高いです！
変数に関して共通因子で説明できる部分が大きいという感じ。

因子負荷量と共通性・独自性を併記するバージョンです。

# 因子負荷量と共通性・独自性の併記バージョン
print('因子負荷量・共通性・独自性:')
pd.concat([loadings, com_uni], axis=1).round(3)

【実行結果】

共通性は変数ごとの因子負荷量の二乗和です。

# 共通性は変数ごとの因子負荷量の二乗和
print('共通性は因子負荷量の二乗和:')
(loadings**2).sum(axis=1).rename('因子負荷量の二乗和').to_frame().round(3)

【実行結果】

７．分散説明率（Prop. Var.）と累積説明率（Cum. Var.）

• 何を測る？

  • 各因子が総分散に占める割合、及び複数因子でどれだけ累積的に説明できるか。

• 値の目安

  • １因子あたりの分散説明率：

    • ≳ 0.10（10%）あれば「まず意味のある因子」

    • ≳ 0.20（20%）あれば「情報量大」

  • 累積説明率（3～5因子まで）：

    • ≳ 0.60（60%）：最低限

    • ≳ 0.70（70%）：標準的

    • ≳ 0.80（80%）：良好

### 7. 因子分散 (因子分散, 分散説明率, 累積説明率) の取得 ---
variance , prop_var, cum_var = fa.get_factor_variance()
print('因子寄与:')
pd.DataFrame({
    '因子分散': variance ,  # 負荷量の二乗和
    '分散説明率': prop_var,
    '累積説明率': cum_var
    }, index=[f'因子{i+1}' for i in range(len(variance ))]
).round(3)

【実行結果】
第３因子までで累積的に $${87.3\%}$$ 説明できています！

８．まとめ（ChatGPTの回答）

• 分析データ自体の評価

  • KMO → 0.6以上を目指し、できれば 0.7 以上

  • Bartlett → p < 0.05 で「因子化 OK」

• 因子数の判断

  • 固有値 → Kaiser基準で λ > 1 を因子数の上限目安に

  • スクリープロット → 「肘」の手前を補助指標に

  • 平行分析 → ランダムではない「意味のある因子」の数

• 推定した因子負荷量の評価

  • 共通性 → 0.5 以上が望ましい

  • 累積説明率 → 70～80%以上でしっかり要約

📕 数式の参考にした資料・記事
とても心強かったです。
ありがとうございます！

（注）PDFファイルを開くリンクです。

テキストの手法とPythonの相性

テキストの因子分析は統計解析ソフト SPSS に依拠しているようです。
一方で Python の主要な因子分析ライブラリは次のとおりです。
・factor_analyzer
・statsmodels.multivariate.factor.Factor
・scikit-learn の sklearn.decomposition.FactorAnalysis

そしてそして…
テキストの因子分析の結果を Python のライブラリでピタリと再現することがとても難しいです（私のスキルでは）。

ですので、「テキストはテキスト」、「PythonはPython」と割り切って進めます。
「似たような解釈ができればＯＫ。些細な違いは気にしない」という気持ちを強くして、おおらかな気持ちで因子分析を楽しみます。

記事の最後はChatGPTが締めくくります。
今回はキャンバスの前で呼吸を整えて。

📘 ChatGPTのひとこと：

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。