「入門はじめての多変量解析」をPythonで写経 Vol.18 ～ ４章「はじめての因子分析」③主因子法

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

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

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

多変量解析の入門を Python と一緒に学ぶ写経シリーズです。
因子分析のうち「探索的因子分析」を取り扱います。

この記事は、因子分析の 因子負荷量の推定 を実践します。
具体的には書籍掲載の 主因子法 の計算を Python の力を借りて取り組みます。

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

はじめに

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

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

４章 はじめての因子分析

この記事は４章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

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

### インポート

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

# 統計処理
import statsmodels.api as sm

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

主因子法による因子分析の結果

先に結論を見ます！
テキストの主因子法によるSPSSの出力結果を Python で追いかけます！
前回記事の振り返り＋主要な出力のご紹介です。

■ 分析の準備
テキスト p.138 表 4.4.1 と同等の表 4.1.2 データをお借りします。

### アンケート調査の結果 p.126 表4.1.2（p.138 表4.4.1と同じ）

# データの登録
data1 = pd.DataFrame(
    {'仕事に集中できる': [3, 4, 3, 1, 2, 5, 1, 4, 2, 5],
     '疲れがとれない': [1, 1, 4, 4, 5, 2, 5, 2, 3, 3],
     'イライラする': [2, 1, 5, 4, 5, 1, 4, 3, 3, 2]},
     index=range(1, 11))
data1.index.name = '被験者No.'

# 結果の表示
data1

【実行結果】
10 名から得た３項目のアンケートの回答データです。

■ データの相関行列
後から使うので計算します。

### データの相関行列
data1.corr()

【実行結果】

■ データの標準化
因子分析のデータは標準化するのが通常のようです。
分母が $${N-1}$$ の標準偏差を用います。

### データの標準化
X_std = (data1 - data1.mean()) / data1.std(ddof=1)

■ 因子分析の実行
statsmodels の Factor() で因子分析を行います。
引数 method='pa' で主因子法を設定します。
実行後、結果の概要（サマリー）を表示します。

### 因子分析の実行 ※テキストp.138 のSPSSの各値とほぼ同じになる

# 因子分析の実行 ※smc=True, maxiter=12を設定
result1_sm = sm.multivariate.Factor(
    X_std, method='pa', smc=True).fit(maxiter=12)

# 結果表示
result1_sm.summary()

【実行結果】
サマリーが表示されました。
最後のテーブル「Pre-rotated loadings」が因子負荷量です！

【引数の補足説明】

• smc
  SMC（※）は「共通性」の下限の推定値です。
  smc=True で共通推定の初期値に SMC を用いる、らしいです。
  ※ ある変数と他の変数との間の重相関係数の二乗

• maxiter
  主因子法の繰り返し回数の上限です。
  テキストとの関連では「12回の反復」の反復回数の上限です。

【出力結果の補足説明】

• 因子分析の固有値：Eigenvalues
  因子負荷量の二乗和であり、因子分散とも呼ばれます。
  変数の分散のうちどれだけ説明しているかの指標になります。

• 共通性：Communality
  変数が共通因子によって説明される程度を表す指標になります。
  因子負荷量の二乗です。

• 因子負荷量（回転前）：Pre-rotated loadings
  お待たせしました！推定したいターゲットです！
  変数が共通因子から受ける影響の度合いを表す指標です。

■ statsmodels から得られる各種指標をご紹介
以下の指標は statamodels の因子分析実行結果から取得できるものたちです。
実際の分析の場で活用されることでしょう。

✅️ 因子負荷量

# 因子負荷量 ※SPSSの「因子行列」に相当
result1_sm.loadings_no_rot

【実行結果】

（参考：テーブル表示）

# 因子負荷量 ※テーブル表示
result1_sm.summary().tables[7].round(3)

【実行結果】

✅️ 共通性

# 共通性 ※SPSSの「共通性－因子抽出後」に相当
result1_sm.communality

【実行結果】

（参考：テーブル表示）

# 共通性 ※テーブル表示
result1_sm.summary().tables[4].round(3)

【実行結果】

✅️ 独自性（１－共通性）

# 独自性: D の対角成分
result1_sm.uniqueness

【実行結果】

（参考：テーブル表示）

# 独自性 ※テーブル表示
1 - result1_sm.summary().tables[4].round(3)

【実行結果】

✅️ 固有値

# 固有値 ※SPSSの「抽出後の負荷量平方和の固有値」に相当
result1_sm.eigenvals

【実行結果】

（参考：テーブル表示）

# 固有値 ※テーブル表示
result1_sm.summary().tables[1].round(3)

【実行結果】

✅️ モデルがフィットした分散共分散行列

# モデルがフィットしたΣ ≒ データの相関行列
result1_sm.fitted_cov

【実行結果】
さきほど計算したデータの相関行列と微妙に異なります。

✅️ 因子得点（Bartlett法）

# 因子得点（Bartlett法） ※statsmodelsのデフォルトはBartlett法
result1_sm.factor_scoring(method='bartlett')

【実行結果】

✅️ 因子得点（回帰法）

# 因子得点（回帰法）
result1_sm.factor_scoring(method='regression')

【実行結果】

（参考：検算）
標準化データ @ 相関行列の逆行列 @ 因子負荷量＝回帰法の因子得点

# 因子得点（回帰法） ※標準化データ @ データの相関行列の逆行列 @ 因子負荷量
(X_std @ np.linalg.inv(X_std.corr()) @ result1_sm.loadings).values

【実行結果】

主因子法の因子負荷量は固有値と固有ベクトルから

主因子法による因子負荷量の推定に進んでいきます！

前回記事の公式から、データの分散共分散行列 $${\Sigma}$$、因子負荷行列 $${\Lambda_f}$$、誤差分散の対角行列 $${D}$$ の関係式を再掲します。

少し式変形します。

分散共分散行列に例題データの相関行列を代入して、その他の項の行列も要素を追加して、テキスト p.139 ～ 140 の式展開を追いかけます。

テキスト p.140 の ☆ の式に追いつきました。

左辺を固有値・固有ベクトルで表現して、テキスト p.141 ～ 143 の式展開を追いかけます（ところどころスキップしてます）。

因子負荷量 $${a_1, a_2, a_3}$$ は「データの相関行列から誤差分散の対角行列を引いた行列」の最大固有値 $${\lambda_1}$$ と固有ベクトル $${\begin{bmatrix}p_1, p_2, p_3\end{bmatrix}^{\top}}$$ との関係で次のようになることが分かりました。

主因子法の計算手続き

因子負荷量の二乗値 $${a_1^2, a_2^2, a_3^2}$$ は、次の「データの相関行列から誤差分散の対角行列を引いた行列」の対角成分 $${1 - \text{Var}(\varepsilon_1), 1 - \text{Var}(\varepsilon_2), 1 - \text{Var}(\varepsilon_3)}$$と一致します。

そこでテキストは、主因子法による因子負荷行列の計算手続きとして、以下の処理を因子負荷量が収束するまで繰り返す方法を p.144 ～ 153 にわたって説明しています。

• 対角成分に「因子負荷量（の途中推定値）の二乗値」を置いてできた行列の固有値・固有ベクトルを算出

• 算出した固有値と固有ベクトルから因子負荷量（の途中推定値）の二乗値を算出して上の行列の対角成分に置く

数式形式で処理の繰り返しイメージを書いてみます。

【行列の更新】
因子負荷量の推定値の二乗を対角成分に設定する

🔽 行列の更新と因子負荷量の推定値の算出を行き来する 🔼

【因子負荷量の推定値の算出】
更新後の行列の固有値・固有ベクトルから因子負荷量の推定値を算出する

繰り返し処理するのはプログラムの得意領域です！
Python で主因子法を実装しましょう。

🔢 Python で主因子法を実装
主因子法による因子負荷量の算出関数を定義します。

### 因子負荷算出関数（主因子法、因子数=1） p.144～153

def exec_factor_analysis_factor1(df, n_iter=100, verbose=False):

    ## 反復0回目の初期処理
    # データの相関行列corr_mtxの算出
    corr_mtx_original = df.corr().values         # データの相関行列
    corr_mtx_updates = corr_mtx_original.copy()  # 反復処理で更新する相関行列
    
    ## 因子負荷算出を反復回数n_iter繰り返し処理
    for i in range(n_iter + 1):

        ## 因子負荷の算出
        # 相関行列の固有値と固有ベクトルの算出
        eig_vals, eig_vecs = np.linalg.eig(corr_mtx_updates)
        # 固有値の降順でindexを取得
        sort_idx = eig_vals.argsort()[::-1]
        # 因子負荷の算出: 因子数=1なので固有値の最大値（１個）を取得して計算
        loadings = eig_vals[sort_idx][0]**(1/2) * eig_vecs[:, sort_idx][:, 0]

        ## 処理経過の表示（Verbose=Trueの場合のみ）
        if verbose:
            print(f'{i}回目の反復')
            print('相関行列:')
            print(corr_mtx_updates)
            print('固有値:', eig_vals[sort_idx][0])
            print('固有ベクトル:', eig_vecs[:, sort_idx][:, 0])
            print('因子負荷:', loadings)
            print('-'*50)
        
        ## 相関行列の対角に因子負荷の二乗値を置き換え
        np.fill_diagonal(corr_mtx_updates, loadings**2)
    
    ## 戻り値: 因子負荷、データの相関行列、固有値、固有ベクトル
    return dict(loadings=loadings, corr_mtx=corr_mtx_original,
                eig_vals=eig_vals[sort_idx], eig_vecs=eig_vecs[:, sort_idx])

【実行結果】なし

主因子法を実行します。
繰り返し回数は引数 n_iter で 12 回を設定します。

# 主因子法による因子負荷の推定（因子数=1）
exec_factor_analysis_factor1(data1, n_iter=12, verbose=True)

【実行結果】
処理１回ごとの途中経過を出力しています。
テキスト p.144 ～ 153 の反復０～反復12 の手動計算結果と同一です。

最終的に確定させた因子負荷量は $${a_1=-0.7319, \ a_2=0.8893, \ a_3=0.9555}$$ です。
テキストの最終結果に一致しています。

🔢 statsmodels で主因子法を実装
statsmodels の因子分析も主因子法に対応しています。
テキスト p.153 の手動計算結果に合致する設定で実行しましょう。
繰り返し処理の最大数を $${13}$$ にします。

### 因子分析の実行 ※テキストp.153 の因子負荷量と同じ結果になる

# 因子分析の実行 ※繰り返し数maxiter=13を設定
result1_sm = sm.multivariate.Factor(
    X_std, method='pa', smc=False).fit(maxiter=13)

# 因子負荷行列の表示
result1_sm.summary().tables[7]

【実行結果】
ピッタリ一致しました！

主因子法による因子負荷量の推定は以上です。
楽しかったですね！

記事の最後は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の教科書です。
よかったらぜひ、お試しくださいませ。

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