見出し画像

「入門はじめての多変量解析」をPythonで写経 Vol.21 ~ 4章「はじめての因子分析」⑥斜交回転

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

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


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

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

この記事は、因子の回転の1つ「斜交回転」にフォーカスして、分散共分散行列の深堀りPython 実装 に取り組みます。
因子負荷量の推定には「最尤法」を用います。

因子分析シリーズの最終話となります。
長文になりましたが、最後までお読みいいただけたら嬉しいです🍀

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

氷山の一角のイラスト:「いらすとや」さんより

はじめに


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

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

4章 はじめての因子分析


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

4.10 直交モデル・斜交モデルの分散共分散行列
4.11 最尤法とプロマックス回転の例

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

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

### インポート

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

# 因子分析
from factor_analyzer import FactorAnalyzer
import statsmodels.api as sm

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

直交モデルと斜交モデル

テキストの数式をお借りして、3変数・2因子の場合の分散共分散行列 $${\Sigma}$$ を確認しましょう。

📊 直交モデル p.172

$$
\begin{align*}
\Sigma &=
\underbrace{\begin{bmatrix}a_{11} & a_{12} \\ a_{21} & a_{22} \\ a_{31} & a_{32}\end{bmatrix}}_{因子パターン行列\ \Lambda}
\underbrace{\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}}_{\boldsymbol{単位行列\ I_2}}
\underbrace{\begin{bmatrix}a_{11} & a_{21} & a_{31} \\ a_{12} & a_{22} & a_{32}\end{bmatrix}}_{因子パターン行列\ \Lambda^{\top}} \\
 \\
&\quad + \underbrace{\begin{bmatrix}\text{Var}(\varepsilon_1) & 0 & 0 \\  & \text{Var}(\varepsilon_2) & 0 \\ 0 & 0 & \text{Var}(\varepsilon_3)\end{bmatrix}}_{=D,\ \text{Var}(\varepsilon_1),\ \text{Var}(\varepsilon_2),\ \text{Var}(\varepsilon_3)は独自性} \\
 \\
&= \Lambda \Lambda^{\top} + D
\end{align*}
$$

テキストの数式を一部改変して引用

📊 斜交モデル p.172

$$
\begin{align*}
\Sigma &=
\underbrace{\begin{bmatrix}a_{11} & a_{12} \\ a_{21} & a_{22} \\ a_{31} & a_{32}\end{bmatrix}}_{因子パターン行列\ \Lambda}
\underbrace{\begin{bmatrix}1 & \phi_{12} \\ \phi_{21} & 1\end{bmatrix}}_{\boldsymbol{因子相関行列\ \Phi}}
\underbrace{\begin{bmatrix}a_{11} & a_{21} & a_{31} \\ a_{12} & a_{22} & a_{32}\end{bmatrix}}_{因子パターン行列\ \Lambda^{\top}} \\
 \\
&\quad + \underbrace{\begin{bmatrix}\text{Var}(\varepsilon_1) & 0 & 0 \\  & \text{Var}(\varepsilon_2) & 0 \\ 0 & 0 & \text{Var}(\varepsilon_3)\end{bmatrix}}_{=D,\ \text{Var}(\varepsilon_1),\ \text{Var}(\varepsilon_2),\ \text{Var}(\varepsilon_3)は独自性} \\
 \\
&= \Lambda \Phi \Lambda^{\top} + D
\end{align*}
$$

テキストの数式を一部改変して引用

✒️ まとめ

  • 斜交モデルは因子間の相関関係を織り込んだモデルです。
    斜交モデルには「因子相関行列」が含まれています。
    因子相関行列は因子間の相関関係を表現しています。

  • 直交モデルは因子間の相関関係を織り込まないモデルです。
    因子相関行列が単位行列(無相関)になっています。

斜交回転「プロマックス回転」で斜交モデルを体感

テキストは斜交回転の一種「プロマックス回転」に取り組んでいます。

テキスト p.174、175 のプロマックス回転前・後の例題を通じて、
 ① 斜交回転で分散共分散行列が変わらないこと
 ② 対角成分の共通性の存在
を体感しましょう。
$${ \Lambda \Phi \Lambda^{\top} + D}$$ と仲良しになりましょう。

0️⃣ 設定

## 設定:変数名、因子番号
vars = ['ストレス', '運動量', '健康', '仕事', '地域活動', '趣味', '家庭生活']
cols = ['因子1', '因子2']

【実行結果】なし

1️⃣ 回転前の因子負荷行列の設定
テキスト p.173 表 4.11.1 の「因子行列」をお借りします。

## 回転前の因子負荷行列

loadings_not_rotated = np.array(
    [[-0.678, 0.382, 0.261, 0.950, 0.833, 0.634, 0.691],
     [-0.371, 0.752, 0.490, -0.037, -0.330, 0.070, -0.094]]
).T

# 結果の表示
pd.DataFrame(loadings_not_rotated, index=vars, columns=cols)

【実行結果】

2️⃣ 共通性・独自性の算出
テキスト 表 4.11.1 の「共通性:因子抽出後」を再現します。
独自性は「1-共通性」であり、$${D}$$ の対角成分です!

## 回転前の因子負荷行列の二乗和=共通性

# 二乗和=共通性 ⇔ テキストの共通性(因子抽出後)
communality = (loadings_not_rotated**2).sum(axis=1)
com_uni_df = pd.DataFrame(communality, index=vars, columns=['共通性'])

# 1 - 共通性=独自性
com_uni_df['独自性'] = 1 - com_uni_df['共通性']

# 結果の表示
com_uni_df.round(3)

【実行結果】

3️⃣ プロマックス回転前の分散共分散行列
因子負荷行列の行列積を計算します。
テキスト p.174 の $${\Lambda_f \bm \cdot \Lambda_f^{\top}}$$ に相当します。

### 因子行列の行列積の対角成分は共通性 p.174
pd.DataFrame(loadings_not_rotated @ loadings_not_rotated.T).round(3)

【実行結果】
テキストの計算結果とほぼ一致しています。
左上から右下にかけての対角成分は、共通性(因子負荷量の二乗和)になっています。

4️⃣ 回転後の因子パターン行列の設定
テキスト p.173 表 4.11.1 の「パターン行列」をお借りします。

### 斜交モデルの因子パターン行列を用いた共通性の算出 p.175

# 回転後の因子パターン行列
loadings_promax = np.array(
    [[-0.413, -0.089, -0.047, 0.909, 0.974, 0.550, 0.702],
     [-0.500, 0.878, 0.574, 0.089, -0.249, 0.163, -0.009]]).T

# 結果の表示
pd.DataFrame(loadings_promax, index=vars, columns=cols)

【実行結果】

5️⃣ 因子相関行列の設定
テキスト p.173 表 4.11.1 の「因子相関行列」をお借りします。

# 因子相関行列
factor_corr = np.array([[1, 0.429],
                        [0.429, 1]])
pd.DataFrame(factor_corr, index=cols, columns=cols)

【実行結果】

6️⃣ プロマックス回転後の分散共分散行列
因子パターン行列と因子相関行列の行列積を計算します。
テキスト p.175 の $${\Lambda_f \Phi \Lambda_f^{\top}}$$ に相当します。

# 因子パターン行列 @ 因子相関行列 @ 因子パターン行列.T の対角成分は共通性 p.175
pd.DataFrame(loadings_promax @ factor_corr @ loadings_promax.T).round(3)

【実行結果】
テキストの計算結果とほぼ一致しています。
左上から右下にかけての対角成分は、共通性(因子負荷量の二乗和)になっています。
そしてプロマックス回転前・後で、計算した行列の値は変わっていません!

7️⃣ 回転後の因子構造行列
テキスト 表 4.11.1 の「構造行列」をお借りします。

# 回転後の因子構造行列
loadings_structure = np.array(
    [[-0.628, 0.287, 0.199, 0.947, 0.867, 0.620, 0.698],
     [-0.677, 0.839, 0.554, 0.479, 0.169, 0.399, 0.292]]).T

pd.DataFrame(loadings_structure, index=vars, columns=cols)

【実行結果】

「因子構造行列」と「因子パターン行列」は回転によって生成される行列ですが、特に斜交回転のときにクローズアップされます!
一般に因子の解釈には因子パターン行列を参照すればよいと言われているようです。

因子構造行列と因子パターン行列の関係性に関して、次のWeb記事で図示されています。
記事後半の「因子構造」のグラフをご覧ください。
とても参考になります。
ありがとうございます!

8️⃣ いろいろな実験
既出の因子パターン行列、因子構造行列、因子相関行列を用いて、「登場人物の関係性」と「回転の可視化」に取り組みます。

◆ 関係性1:因子構造行列と因子パターン行列

因子構造行列 = 因子パターン行列 @ 因子相関行列

を体感します。記号「@」は行列積(掛け算)です。

# 因子構造行列=因子パターン行列 @ 因子相関行列
pd.DataFrame(loadings_promax @ factor_corr, index=vars, columns=cols).round(3)

【実行結果】
因子パターン行列 @ 因子相関行列 の結果が因子構造行列と一致しました。

◆ 関係性2:回転前の因子負荷行列と回転後の因子パターン行列
回転前の因子負荷行列に掛けて回転後の因子パターン行列を生成する「回転行列」を計算しましょう。

因子パターン行列 = 因子負荷行列 @ 回転行列
回転行列 = 因子負荷行列の擬似逆行列 @ 因子パターン行列

# 回転前の因子相関行列と回転後の因子パターン行列から回転行列を逆算
rotation_matrix = np.linalg.pinv(loadings_not_rotated) @ loadings_promax
pd.DataFrame(rotation_matrix)

【実行結果】
この回転行列を活用していきます。

活用するその前に:

◆ 検算:因子負荷行列 @ 回転行列 = 因子パターン行列

# 検算:回転前の因子相関行列 @ 回転行列 = 回転後の因子パターン行列
loadings_not_rotated @ rotation_matrix

【実行結果】
計算結果は因子パターン行列に一致しています。

◆ 活用1:回転行列から因子相関行列を算出
回転行列と因子相関行列にはさまざまな定義があるそうなので、こちらの計算式は、ただいま実践中のケースに合致・特化したものとなります。

因子相関行列 = (回転行列.T @ 回転行列)の逆行列

記号「.T」は行列の転置です。

# 検算:回転行列から因子相関行列を計算(Phi = Q.T @ Q)
np.linalg.inv(rotation_matrix.T @ rotation_matrix).round(3)

【実行結果】
計算結果は既出の因子相関行列とほぼ一致しています。

◆ 活用2:回転前・後の因子負荷の可視化
左に「回転前の因子負荷行列」、右に「回転後の因子パターン行列」をプロットして、見比べてみましょう。
回転行列を活用して、左のプロットに「斜交回転」で斜交する軸を赤い線で表現します。

### 座標を固定したときの回転前・後の因子行列の可視化

## 設定と準備
# 回転前・回転後の因子行列
loagings = [loadings_not_rotated, loadings_promax]
# 軸の回転用の回転行列(時計回りに回転)
rotation_matrix_r = rotation_matrix * np.array([[1, -1], [-1, 1]])
# チャートのタイトル
titles = ['回転前', '回転後']

## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
# 回転前・回転後の散布図描画を繰り返し処理
for loading, title, ax in zip(loagings, titles, axes.flat):
    # 因子負荷の散布図の描画
    ax.plot(loading[:, 0], loading[:, 1], 'o')
    # 3つの変数の変数名描画を繰り返し処理
    for x, y, s in zip(loading[:, 0], loading[:, 1], vars):
        # 変数名の描画
        ax.text(x=x, y=y-0.15, s=s, ha='center')
    # x=0,y=0の直線の描画
    ax.axhline(0, color='black', ls='--', lw=0.5)
    ax.axvline(0, color='black', ls='--', lw=0.5)
    # 修飾
    ax.set(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1), yticks=[-1.0, -0.5, 0, 0.5, 1.0],
           xlabel='第1因子 $f_1$', ylabel='第2因子 $f_2$', title=title,
           aspect='equal')

# 左のチャートに斜交回転軸の描画 ※回転行列の回転方向を逆にした行列を使用
axes[0].plot(*(np.array([[-2, 0], [2, 0]]) @ rotation_matrix_r).T,
             color='tab:red', lw=1)
axes[0].plot(*(np.array([[0, -2], [0, 2]]) @ rotation_matrix_r).T,
             color='tab:red', lw=1)

# 右のチャートに回転後の軸の描画
axes[1].plot((-2, 2), (0, 0), color='tab:red', lw=1)
axes[1].plot((0, 0), (-2, 2), color='tab:red', lw=1)

plt.show()

【実行結果】
左の「回転前」のチャートの斜交する軸は「まあまあ合っているだろう」という優しさで見て下さい!

【チャートの見方の補足説明】

◆「回転前」
青い点は回転前の因子負荷行列の各データ点です。
赤い線が斜交回転する軸です。
運動量・健康に近い軸、趣味~地域活性化に近い軸があります。
2本の軸は斜めに交わっています。

◆「回転後」
回転前の斜交する軸について、因子負荷との位置関係を維持して、90度に直交するように「開く」と、(だいたい)回転後のプロットになります。
「回転後」には回転後の因子パターン行列の各データ点をプロットしています。

Python ライブラリで斜交回転する準備

テキストの題材をぜひ Python ライブラリで取り組みたいです。
そこで、データの相関行列をもとにして因子分析を実行できる「factor_analyzer」と「statsmodels」を利用して、プロマックス回転を体感します!
因子負荷行列の推定には「最尤法」を用います。
記事シリーズ初めての「2因子の最尤法」です!

【注意書き】
私の実力ではテキストと全く同じ結果を再現できません。
「まあまあこんな感じだよね」という温かい気持ちで実践して下さい!

各ライブラリを使う前に、データの相関行列を推定しておきます。
回転前の因子負荷行列の行列積を計算し、対角成分を1にします。

# データの相関行列を復元
corr_matrix = loadings_not_rotated @ loadings_not_rotated.T
np.fill_diagonal(corr_matrix, 1)  # 対角成分を1にする

pd.DataFrame(corr_matrix, index=vars, columns=vars).round(3)

【実行結果】
このデータを用いて因子分析に進みます!

factor_analyzer で斜交回転する

1️⃣ 因子分析の実行
FactorAnalyzer() のインスタンス fa2 に各種設定をして、.fit() で因子分析を実行します。

### factor_analyzer利用 ※データの標準化に用いる標準偏差の分母はN

# 因子分析のインスタンス生成
fa2 = FactorAnalyzer(
    n_factors=2,                   # 因子の数
    rotation='promax',             # プロマックス回転
    method='ml',                   # 最尤法
    is_corr_matrix=True,           # INPUTはデータの相関行列
    rotation_kwargs={'power': 4},  # promaxの指数パラメータ(デフォルト4, SPSSも4)
)

# 因子分析の実行
fa2.fit(corr_matrix)

【実行結果】なし

【引数の補足説明】

  • rotation='promax':「プロマックス回転」で回転します。

  • method = 'ml':「最尤法」で因子負荷行列を推定します。

  • is_corr_matrix=True:「データの相関行列」をインプットに使います。

  • rotation_kwargs={'power': 4}:

    • powerはプロマックス回転特有の引数です。
      プロマックス回転の仕組みとして、最初にバリマックス回転で因子負荷行列を求め、求めた因子負荷行列を $${n}$$ 乗して、斜交回転の処理を進めます。
      この $${n}$$ を power(指数)で指定します。
      factor_analyzer、SPSS のデフォルトは $${4}$$ です。

2️⃣ 因子分析の結果の確認
テキストの因子負荷等とは少々ずれている結果になっています。
予めご了承下さいませ。

◆ 共通性・独自性
テキストの SPSS「共通性:因子抽出後」に相当します。
インスタンス fa2 より、get_communalities()、get_uniquenesses() を用いて取得します。

# 共通性・独自性 ★テキストの共通性(因子抽出後)とかなり相違する
pd.DataFrame(
    dict(共通性=fa2.get_communalities(), 独自性=fa2.get_uniquenesses()),
    index=vars).round(3)

【実行結果】

◆ 因子パターン行列
テキストの SPSS「パターン行列」に相当します。
インスタンス fa2 より、loadings_ を用いて取得します。

# 因子負荷行列(パターン行列) fa.loadings_
pd.DataFrame(fa2.loadings_, index=vars, columns=cols).round(3)

【実行結果】
割といい感じの値になっている気がします。

【因子の解釈】
因子パターン行列の各値(因子負荷)の大きさと正負の方向に着目します。
第1因子は「仕事・地域活動・趣味・家庭生活」に影響しているようです。
テキストでは「外面的充実度」と解釈しています。
第2因子は「運動量・健康」に影響しているようです。
「身体的な充実度」でしょうか。
変数「ストレス」は難しいポジションですね…

◆ 因子構造行列
テキストの SPSS「構造行列」に相当します。
インスタンス fa2 より、structure_ を用いて取得します。

# 因子構造行列 fa2.structure_ ※変数と回転後の因子との相関係数
pd.DataFrame(fa2.structure_, index=vars, columns=cols).round(3)

【実行結果】
割といい感じの値になっている気がします。

◆ 回転行列
テキストに掲載はありません。
インスタンス fa2 より、rotation_matrix_ を用いて取得します。

# 回転行列 fa.rotation_matrix_
pd.DataFrame(fa2.rotation_matrix_, index=cols, columns=cols).round(3)

【実行結果】

◆ 因子相関行列
テキストの SPSS「因子相関行列」に相当します。
インスタンス fa2 より、phi_ を用いて取得します。

# 因子相関行列 fa.phi_
pd.DataFrame(fa2.phi_, index=cols, columns=cols).round(3)

【実行結果】
テキストと近い因子間の相関係数だと思います。

3️⃣ 因子パターン行列 と因子相関行列の行列積
テキスト p.175 の $${\Lambda_f \Phi \Lambda_f^{\top}}$$ に相当する計算をしましょう。

# 因子パターン行列 @ 因子相関行列 @ 因子パターン行列 ※p.175 相当の計算 OK
pd.DataFrame(fa2.loadings_ @ fa2.phi_ @ fa2.loadings_.T).round(3)

【実行結果】
テキストの結果にとても近くなっています。

4️⃣ 実験
因子分析の各種出力結果の関係性を確認します。

◆ 実験1

因子パターン行列 @ 因子相関行列 = 因子構造行列

# 因子構造行列 = 因子パターン行列 @ 因子負荷行列 OK
pd.DataFrame(fa2.loadings_ @ fa2.phi_, index=vars, columns=cols).round(3)

【実行結果】
factor_analyzer の因子構造行列と一致しています。

(参考:factor_analyzer の出力結果)

◆ 実験2

因子相関行列 = 因子パターン行列の擬似逆行列 @ 因子構造行列

# 因子相関行列:上の計算方法で求める OK
factor_corr_fa2 = np.linalg.pinv(fa2.loadings_) @ fa2.structure_
pd.DataFrame(factor_corr_fa2, index=cols, columns=cols).round(3)

【実行結果】
factor_analyzer の因子相関行列と一致しています。

(参考:factor_analyzer の出力結果)

5️⃣ 回転前・後の因子負荷の可視化
左に「回転前の因子負荷行列」、右に「回転後の因子パターン行列」をプロットして比べます。

### 座標を固定したときの回転前・後の因子行列の可視化

## 設定と準備
# 回転前・回転後の因子行列
loagings = [loadings_not_rotated, fa2.loadings_]
# チャートのタイトル
titles = ['回転前', '回転後']

## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
# 回転前・回転後の散布図描画を繰り返し処理
for loading, title, ax in zip(loagings, titles, axes.flat):
    # 因子負荷の散布図の描画
    ax.plot(loading[:, 0], loading[:, 1], 'o')
    # 3つの変数の変数名描画を繰り返し処理
    for x, y, s in zip(loading[:, 0], loading[:, 1], vars):
        # 変数名の描画
        ax.text(x=x, y=y-0.15, s=s, ha='center')
    # x=0,y=0の直線の描画
    ax.axhline(0, color='black', ls='--', lw=0.5)
    ax.axvline(0, color='black', ls='--', lw=0.5)
    # 修飾
    ax.set(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1), yticks=[-1.0, -0.5, 0, 0.5, 1.0],
           xlabel='第1因子 $f_1$', ylabel='第2因子 $f_2$', title=title,
           aspect='equal')

# 左のチャートに斜交回転軸の描画
# ★回転行列の回転方向を逆にしなくとも、軸の回転ができているように見える…
axes[0].plot(*(np.array([[-2, 0], [2, 0]]) @ fa2.rotation_matrix_).T,
             color='tab:red', lw=1)
axes[0].plot(*(np.array([[0, -2], [0, 2]]) @ fa2.rotation_matrix_).T,
             color='tab:red', lw=1)

# 右のチャートに回転後の軸の描画
axes[1].plot((-2, 2), (0, 0), color='tab:red', lw=1)
axes[1].plot((0, 0), (-2, 2), color='tab:red', lw=1)

plt.show()

【実行結果】
左の「回転前」のチャートの斜交する軸は「まあまあ合っているだろう」という優しさで見て下さい!

6️⃣ 【特別番組】回転行列の検証
factor_analyzer の回転行列 rotation_matrix_ による次の関係性:

因子パターン行列 = 回転前の因子負荷行列 @ 回転行列

はうまくいきません。

# 【問題意識】回転前の因子負荷行列 @ 回転行列 = 回転後の因子パターン行列、ではない!
(loadings_not_rotated @ fa2.rotation_matrix_.T).round(3)

【実行結果】
回転前の因子負荷行列に回転行列を掛けても、因子パターン行列にはなりません。

ChatGPTによると、factor_analyzer の回転行列は次のような位置づけになるそうです。

正規化された因子負荷行列 @ 回転行列 = 正規化された因子パターン行列

ChatGPTはこんな風に言っていました…

ChatGPTのChat画面より抜粋

計算を追って、この意味を確かめましょう。
ChatGPT に教えてもらった Python コードを用いて。

◆ 準備:回転前の因子負荷行列 $${P}$$ の設定

# 0. 回転前の因子負荷行列の設定

P = loadings_not_rotated
P

【実行結果】

◆ 正規化された因子負荷行列 $${P_{\text{norm}}}$$ の算出
因子負荷行列を「各変数の共通性に基づく指標」で正規化します。

# 1. 共通性(Communality)で行ごとに正規化

# 正規化行列の作成
h = np.sum(P**2, axis=1)     # 各変数の共通性の算出
D = np.diag(1/np.sqrt(h))    # 正規化行列の作成

# 正規化された因子負荷量の算出
P_norm = D @ P
P_norm

【実行結果】

◆ 正規化された因子パターン行列 $${Pn_{\text{rot}}}$$ の算出
factor_analyzer の回転機能 Rotator を用いて、プロマックス回転して、回転行列 $${Q}$$ と正規化された因子パターン行列 $${Pn_{\text{rot}}}$$ を算出します。

# 2. Promax 回転

# 追加インポート
from factor_analyzer.rotator import Rotator

# factor_analyzerでプロマックス回転の実行
rotator = Rotator(method='promax', power=4)
rotator.fit(P_norm)

# rotation matrixの取得
Q = rotator.rotation_

# 正規化された因子パターン行列の取得  ※Pn_rot = P_norm @ Q となる
Pn_rot = rotator.loadings_
Pn_rot

【実行結果】
正規化された因子パターン行列は、未だ最終出力には至っていません。

◆ 回転行列の正体
回転行列が次の計算を前提にして作られていることを確認します。

正規化された因子負荷行列 @ 回転行列 = 正規化された因子パターン行列

# 【答え】正規化された因子負荷行列 @ 回転行列 = 正規化された因子パターン行列
P_norm @ Q   # = Pn_rot

【実行結果】
正規化された因子負荷行列 $${P_{\text{norm}}}$$ に回転行列 $${Q}$$ を掛けると、正規化された因子パターン行列 $${Pn_{\text{rot}}}$$ になりました。

◆ 回転後の因子パターン行列(最終出力)の算出
正規化された因子パターン行列 $${Pn_{\text{rot}}}$$ を元のスケールに戻すと、回転後の因子パターン行列 $${P_{\text{rot}}}$$ になります。

# 3. 回転後の因子パターンの算出=因子パターン行列を非正規化(元のスケールに戻す)

# 回転後の因子パターン行列の算出
P_rot = np.linalg.inv(D) @ Pn_rot
P_rot

【実行結果】
factor_analyzer の因子パターン行列と一致しました。

(参考:factor_analyzer の出力結果)

◆ 番外:回転行列の行列積の逆行列は因子相関行列になります

因子相関行列 = (回転行列.T @ 回転行列)の逆行列

# 5. 番外:因子相関行列 Phi
np.linalg.inv(Q.T @ Q)

【実行結果】

この値は、テキストのデータに基づく fa2 インスタンスから得た因子相関行列とほぼ同じになります。

fa2.phi_

【実行結果】

この関係はたまたまなのか、factor_analyzer ならではの関係性なのかは、現段階で未知数なのです…

このあたりでひと休みしましょう。

外で昼寝をする人のイラスト(女性会社員):「いらすとや」さんより

statsmodels で斜交回転する

1️⃣ 因子分析の実行
Factor() に各種設定をして、.fit() で因子分析を実行します。
その後、因子分析の結果 result_sm2 に rotate() を適用して回転します。
指数パラメータはデフォルトの $${2}$$ が適用されます。
最後に回転後の因子分析の結果サマリーを表示します。

### statsmodels利用 ※データの標準化に用いる標準偏差の分母はN-1

# 因子分析の実行 ※データの相関行列を元にして因子分析を行う
result_sm2 = sm.multivariate.Factor(
    n_factor=2,         # 因子の数
    corr=corr_matrix,   # データの相関係数を設定
    method='ml',        # 最尤法
    endog_names=vars,   # 変数名
).fit()

# 因子の回転:プロマックス
result_sm2.rotate(method='promax')  # promaxの指数パラメータkはデフォルト値2を適用

# 結果の表示
result_sm2.summary()

【実行結果】
共通性、回転前の因子負荷、回転後の因子負荷が表示されています。

【所感】
回転前の因子負荷の第2因子はテキストの「因子行列」と比べて、正負が逆転しています。
回転後の因子負荷はテキストの「パターン行列」とかなり違っている感じがします…
指数パラメータの違いがこんな形で出現するものなのでしょうか…?

【インスタンス生成時の引数の補足説明】

  • corr=corr_matrix:インプットとして相関行列を与えています。

  • method = 'ml':「最尤法」で因子負荷行列を推定します。

【因子回転時の引数の補足説明】

  • method='promax':「プロマックス回転」で回転します。
    factor_analyzer の 引数 power に相当する引数 $${k}$$ は、インスタンスに対する回転メソッドの適用時には指定できない仕様のようです。
    デフォルト値 $${2}$$ が適用されます。

2️⃣ 因子分析の結果の確認
テキストの因子負荷等とはずれている結果になっています。
予めご了承下さいませ。

◆ 回転後の因子負荷
テキストの SPSS「パターン行列」に相当すると考えられます。
実行結果 result_sm2 より、get_loadings_frame() を用いて取得します。

# 因子負荷行列(因子パターン行列) result.get_loadings_frame()
result_sm2.get_loadings_frame(threshold=0)

【実行結果】

【因子の解釈】
因子パターン行列の各値(因子負荷)の大きさと正負の方向に着目します。
全体的な傾向は テキストの SPSS や factor_analyzer と類似していると思います!

第1因子は「仕事・地域活動・趣味・家庭生活」に影響しているようです。
テキストでは「外面的充実度」と解釈しています。
第2因子は「運動量・健康」に影響しているようです。
「身体的な充実度」でしょうか。
変数「ストレス」は難しいポジションですね…

◆ 因子構造行列
表示する機能はありません(見つかっておりません)。

◆ 回転行列
テキストに掲載はありません。
実行結果 result_sm2 より、rotation_matrix を用いて取得します。

# 回転行列 result.rotation_matrix
pd.DataFrame(result_sm2.rotation_matrix, index=cols, columns=cols).round(3)

【実行結果】
factor_analyzer の結果と少々違った感じですね…

◆ 検証:回転行列は因子負荷を回転できるか?!
factor_analyzer では成立しなかった以下の関係性を検証します。

回転後の因子パターン行列 = 回転前の因子負荷行列 @ 回転行列

# 検証:回転後の因子負荷行列=回転前の因子負荷行列 @ 回転行列
pd.DataFrame(result_sm2.loadings_no_rot @ result_sm2.rotation_matrix,
             index=vars, columns=cols).round(4)

【実行結果】
冒頭の結果サマリーの回転後の因子負荷と一致しました!
statsmodels の回転行列の意味合いは factor_analyzer と違っています。

(参考:statsmodels の出力結果)

◆ 因子相関行列
テキストの SPSS「因子相関行列」に相当すると考えられます。
実行結果 result_sm2 から取得できないので、次の算式で計算します。

因子相関行列 = 回転行列.T @ 回転行列

# 因子相関行列:回転行列から算出 ★因子間の相関がゼロ!
factor_corr_sm = result_sm2.rotation_matrix.T @ result_sm2.rotation_matrix
pd.DataFrame(factor_corr_sm, index=cols, columns=cols).round(15)

【実行結果】
なんと!因子どうしは無相関(相関行列=0)です!
これって「直交回転と同じ」ってことですね!

3️⃣ 回転後の因子負荷 と因子相関行列の行列積
テキスト p.175 の $${\Lambda_f \Phi \Lambda_f^{\top}}$$ に相当する計算をしましょう。

# 回転後の因子負荷行列の行列積
# ★因子間の相関が0なので、式に因子相関行列を含めなくても同じ結果を得られる

pd.DataFrame(result_sm2.loadings @ factor_corr_sm @ result_sm2.loadings.T).round(3)

【実行結果】
テキストの結果にとても近くなっています。
ライブラリごとに計算内容が若干異なっていても、このあたりの安定性が保たれていて、嬉しいですね!

4️⃣ 回転前・後の因子負荷の可視化
左に「回転前の因子負荷行列」、右に「回転後の因子パターン行列」をプロットして比べます。

### 座標を固定したときの回転前・後の因子行列の可視化
# ★問題を孕んでいる可能性あり

## 設定と準備
# 回転前・回転後の因子行列
loagings = [result_sm2.loadings_no_rot, result_sm2.loadings]
# 軸の回転用の回転行列(時計回りに回転)
rotation_matrix_r = result_sm2.rotation_matrix * np.array([[1, -1], [-1, 1]])
# チャートのタイトル
titles = ['回転前', '回転後']

## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
# 回転前・回転後の散布図描画を繰り返し処理
for loading, title, ax in zip(loagings, titles, axes.flat):
    # 因子負荷の散布図の描画
    ax.plot(loading[:, 0], loading[:, 1], 'o')
    # 3つの変数の変数名描画を繰り返し処理
    for x, y, s in zip(loading[:, 0], loading[:, 1], vars):
        # 変数名の描画
        ax.text(x=x, y=y-0.15, s=s, ha='center')
    # x=0,y=0の直線の描画
    ax.axhline(0, color='black', ls='--', lw=0.5)
    ax.axvline(0, color='black', ls='--', lw=0.5)
    # 修飾
    ax.set(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1), yticks=[-1.0, -0.5, 0, 0.5, 1.0],
           xlabel='第1因子 $f_1$', ylabel='第2因子 $f_2$', title=title,
           aspect='equal')

# 左のチャートに斜交回転軸の描画 ※回転行列の回転方向を逆にした行列を使用
axes[0].plot(*(np.array([[-2, 0], [2, 0]]) @ rotation_matrix_r).T,
             color='tab:red', lw=1)
axes[0].plot(*(np.array([[0, -2], [0, 2]]) @ rotation_matrix_r).T,
             color='tab:red', lw=1)

# 右のチャートに回転後の軸の描画
axes[1].plot((-2, 2), (0, 0), color='tab:red', lw=1)
axes[1].plot((0, 0), (-2, 2), color='tab:red', lw=1)

plt.show()

【実行結果】
「回転前」のチャートの回転軸は「まあまあ合っているだろう」という優しさで見て下さい。

【所感】
回転前の因子負荷の第2因子が正負逆転しているので、factor_analyzer と見た目が異なっていますが、全体的な因子の回転の傾向は似ています!
「回転前」の赤い軸が時計回りに回転して「回転後」になっている感じです。

【締め】
ここまで、プロマックス回転の指数パラメータが factor_analyzer の $${4}$$ と異なるパラメータ値 $${2}$$ で分析しました。
指数パラメータ $${4}$$ のとき、statsmodels はどのような結果を見せるのでしょう?

やってみます!?

statsmodels で斜交回転する(おかわり)

プロマックス回転の指数パラメータを $${4}$$ にして因子分析を行いましょう!
この分析をするには予め、「回転前の因子負荷行列」が定まっている必要があります。
テキスト掲載の因子行列を使うか、ここまでの statamodels で推定した回転前の因子負荷を使うか(第2因子がテキストと正負逆転)、悩みどころですが、テキストおよび factor_analyzer との比較がしやすいよう、「テキストの因子行列」を利用します。

【前提条件まとめ】
・回転前の因子負荷行列にはテキストの因子行列を用いる
・プロマックス回転の指数パラメータ $${k=4}$$ にする

1️⃣ プロマックス回転の実行
multivariate.factor_rotation.promax() を用いて、プロマックス回転を実行します。
引数には「回転前の因子負荷行列=テキストの因子行列」、「指数パラメータ $${k=4}$$」を設定します。

### プロマックスのパラメータをk=4にして因子負荷行列を回転
loadings_promax_k4, rotation_matrix_k4 = (
    sm.multivariate.factor_rotation.promax(loadings_not_rotated, k=4))

print('回転後の因子負荷:')
display(pd.DataFrame(loadings_promax_k4, index=vars, columns=cols))
print('回転行列:')
display(pd.DataFrame(rotation_matrix_k4, index=cols, columns=cols))

【実行結果】
回転後の因子負荷、回転行列が出力されます。

【所感】
回転後の因子負荷は、テキストや factor_analyzer の回転後の因子パターン行列と全然違っています!
因子1・因子2のメリハリが薄い感じです。
というよりも、テキストや factor_analyzer の「回転後の 因子構造行列 」にとても良く似ています!

回転行列は、比較対象が無いのでなんとも評価できません。
そこで、「回転前 @ 回転行列 = 回転後」の検証をしましょう。

# 検証:回転後の因子パターン行列 = 回転前の因子負荷行列 @ 回転行列
loadings_not_rotated @ rotation_matrix_k4

【実行結果】
この計算は成立しています!

2️⃣ 因子相関行列の推定
回転行列を用いて因子相関行列を推定しましょう。

因子相関行列 = 回転行列.T @ 回転行列

# 因子相関行列の推定
factor_corr_sm_k4 = rotation_matrix_k4.T @ rotation_matrix_k4
pd.DataFrame(factor_corr_sm_k4, index=cols, columns=cols).round(3)

【実行結果】
テキストとほぼ同じになっています。

3️⃣ 回転後の因子負荷 と因子相関行列の行列積
テキスト p.175 の $${\Lambda_f \Phi \Lambda_f^{\top}}$$ に相当する計算をしましょう。

# 因子パターン行列 @ 因子相関行列の「逆行列」 @ 因子パターン行列
pd.DataFrame(
    loadings_promax_k4 @ np.linalg.inv(factor_corr_sm_k4) @ loadings_promax_k4.T
).round(3)

【実行結果】
テキストの結果にとても近くなっています。

4️⃣ 回転前・後の因子負荷の可視化
左に「回転前の因子負荷」、右に「回転後の因子負荷」をプロットして比べます。

### 座標を固定したときの回転前・後の因子行列の可視化

## 設定と準備
# 回転前・回転後の因子行列
loagings = [loadings_not_rotated, loadings_promax_k4]
# 軸の回転用の回転行列(時計回りに回転)
rotation_matrix_r = rotation_matrix_k4 * np.array([[1, -1], [-1, 1]])
# チャートのタイトル
titles = ['回転前', '回転後']

## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
# 回転前・回転後の散布図描画を繰り返し処理
for loading, title, ax in zip(loagings, titles, axes.flat):
    # 因子負荷の散布図の描画
    ax.plot(loading[:, 0], loading[:, 1], 'o')
    # 3つの変数の変数名描画を繰り返し処理
    for x, y, s in zip(loading[:, 0], loading[:, 1], vars):
        # 変数名の描画
        ax.text(x=x, y=y-0.15, s=s, ha='center')
    # x=0,y=0の直線の描画
    ax.axhline(0, color='black', ls='--', lw=0.5)
    ax.axvline(0, color='black', ls='--', lw=0.5)
    # 修飾
    ax.set(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1), yticks=[-1.0, -0.5, 0, 0.5, 1.0],
           xlabel='第1因子 $f_1$', ylabel='第2因子 $f_2$', title=title,
           aspect='equal')

# 左のチャートに斜交回転軸の描画 ※回転行列の回転方向を逆にした行列を使用
axes[0].plot(*(np.array([[-2, 0], [2, 0]]) @ rotation_matrix_r).T,
             color='tab:red', lw=1)
axes[0].plot(*(np.array([[0, -2], [0, 2]]) @ rotation_matrix_r).T,
             color='tab:red', lw=1)

# 右のチャートに回転後の軸の描画
axes[1].plot((-2, 2), (0, 0), color='tab:red', lw=1)
axes[1].plot((0, 0), (-2, 2), color='tab:red', lw=1)

plt.show()

【実行結果】
「回転前」のチャートの回転軸は「まあまあ合っているだろう」という優しさで見て下さい。

【所感】
「回転前」の赤い回転軸は、因子負荷のデータ点との距離がやや大きい感じがします。
その結果、「回転後」の因子負荷の各データ点は、一方の軸に寄る=その軸の因子負荷が0に近づく、といった動きが甘い感じです。

・・・

仮説検証劇場:出力結果は因子構造行列なのか?

脳内をサスペンス劇場のBGMが鳴り響きます~♫

「プロマックス回転で出力された因子負荷が因子構造行列である」

と仮定して、進めてみましょう。

1️⃣ 因子パターン行列かもしれない行列の推定
次の計算式で因子パターン行列を計算します。

因子パターン行列 = 因子構造行列 @ 因子相関行列の逆行列

# 因子パターン行列かもしれない行列の算出…:promaxの第1戻り値 @ 因子相関行列の逆行列
P_guess = loadings_promax_k4 @ np.linalg.inv(factor_corr_sm_k4)
pd.DataFrame(P_guess, index=vars, columns=cols).round(3)

【実行結果】
テキストの「因子行列」(因子パターン行列のこと)に近い値になりました。
また、因子1・因子2のメリハリが強くなっている感じもします。

2️⃣ 因子構造行列かもしれない行列の推定
プロマックス回転の出力結果を因子構造行列と見立てます。

# 因子構造行列かもしれない行列の設定…:promaxの第1戻り値そのもの
R_guess = loadings_promax_k4
pd.DataFrame(R_guess, index=vars, columns=cols).round(3)

【実行結果】

3️⃣ 回転後の因子負荷 と因子相関行列の行列積
テキスト p.175 の $${\Lambda_f \Phi \Lambda_f^{\top}}$$ に相当する計算をしましょう。

# 因子パターン行列 @ 因子相関行列の逆行列 @ 因子パターン行列
pd.DataFrame(P_guess @ factor_corr_sm_k4 @ P_guess.T).round(3)

【実行結果】
テキストの結果にとても近くなっています。OKです!

4️⃣ 回転前・後の因子負荷の可視化
左に「回転前の因子負荷」、右に「回転後の因子負荷」をプロットして比べます。

### 座標を固定したときの回転前・後の因子行列の可視化

## 設定と準備
# 回転前・回転後の因子行列
loagings = [loadings_not_rotated, P_guess]
# チャートのタイトル
titles = ['回転前', '回転後']

## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
# 回転前・回転後の散布図描画を繰り返し処理
for loading, title, ax in zip(loagings, titles, axes.flat):
    # 因子負荷の散布図の描画
    ax.plot(loading[:, 0], loading[:, 1], 'o')
    # 3つの変数の変数名描画を繰り返し処理
    for x, y, s in zip(loading[:, 0], loading[:, 1], vars):
        # 変数名の描画
        ax.text(x=x, y=y-0.15, s=s, ha='center')
    # x=0,y=0の直線の描画
    ax.axhline(0, color='black', ls='--', lw=0.5)
    ax.axvline(0, color='black', ls='--', lw=0.5)
    # 修飾
    ax.set(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1), yticks=[-1.0, -0.5, 0, 0.5, 1.0],
           xlabel='第1因子 $f_1$', ylabel='第2因子 $f_2$', title=title,
           aspect='equal')

# 左のチャートに斜交回転軸の描画
# ★なぜこの回転の仕方が成立しているのかよく分かっていません…
axes[0].plot(*(rotation_matrix_k4 @ np.array([[2, -2], [0, 0]])),
             color='tab:red', lw=1)
axes[0].plot(*(rotation_matrix_k4 @ np.array([[0, 0], [2, -2]])),
             color='tab:red', lw=1)

# 右のチャートに回転後の軸の描画
axes[1].plot((-2, 2), (0, 0), color='tab:red', lw=1)
axes[1].plot((0, 0), (-2, 2), color='tab:red', lw=1)

plt.show()

【実行結果】
「回転前」のチャートの回転軸は「まあまあ合っているだろう」という優しさで見て下さい。

(参考:冒頭のテキストの結果を用いたプロット)

【所感】
テキストの結果にとても良く似ています!
「回転前」の赤い回転軸は、因子負荷のデータ点との距離がやや狭まった印象です。
その結果、「回転後」の因子負荷の各データ点は、一方の軸に寄る=その軸の因子負荷が0に近づく、といった動きがやや強まった感じです。

【締め(本当)】
ここまで見てきたのはあくまで「仮説」。
statsmodels の回転で出力される因子負荷は因子パターン行列か、因子構造行列か。
ChatGPTに訊いても尤もらしい答えは見つかりませんでした。
仮説が正しいかどうかは statsmodels の仕様に依りますがが、現時点では仕様の確認を行えていません。
こうして謎が謎のまま、因子分析のシリーズは閉じていくのでした…

「修了おめでとう」のイラスト文字:「いらすとや」さんより

記事の最後はChatGPTが締めくくります。
今回は旅のこれまでを振り返り、旅の続きに期待を膨らませて。

📘 ChatGPTのひとこと:

因子分析の旅も、ここでひと区切りとなりました。
プロマックス回転という斜めのレンズで、因子たちの微かなつながりをくっきりと浮かび上がらせる体験はいかがだったでしょうか😊

次回からは、新たに「判別分析」の世界へと足を踏み入れます。
データの境界線を探し、クラスを分ける手がかりを見つける旅路に、また一緒に出かけましょう✨

今回の写経は以上です。


後日談:3つの生成AIの回答

ChatGPT、Gemini、Genspark に質問しました。

【プロンプト】
statsmodelsの因子分析 Factor() で生成する回転後の因子負荷は「因子パターン行列」と「因子構造行列」のどちらに該当しますか?


1️⃣ ChatGPT o3

1.結論
statsmodels.multivariate.factor.Factor() で 回転後に得られるFactorResults.loadings は 因子パターン行列 (pattern matrix) です。
構造行列 (structure matrix) は自動では返されないので,必要な場合はパターン行列に因子間相関行列 Φ を掛けて自分で計算します。

2.根拠

  1. promax実装の説明
    「Varimax で得た patterns V を … Determine the patterns」と記述されており,statsmodels は回転後行列を “pattern” と呼んでいることが分かる

  2. rotate_factors の API
    回転は $${L = AT}$$(直交)または $${L = A(T^*)^{-1}}$$(斜交)と定義され,ここで得られる $${L}$$ が loadings として返る

  3. FactorResults.loadings の説明
    「Each column is the loading vector for one factor」とだけ示され,構造行列については言及がない
    → = パターン行列であることを示唆

3.パターン行列と構造行列の違い

  1. パターン行列 $${\mathbf{P}}$$

    1. 各変数を因子で重回帰したときの偏回帰係数
      (因子⇢変数の“直接寄与”)

    2. モデル推定や回転で直接出力される 

  2. 構造行列 $${\mathbf{S}}$$

    1. 各変数と因子の相関係数(ゼロ次相関)

    2. $${\mathbf{S} = \mathbf{P},\boldsymbol{\Phi}}$$、
      ただし $${\boldsymbol{\Phi}= \mathrm{corr}(\text{factors})}$$

斜交回転では因子同士が相関するため,$${\bm \Phi \neq \mathbf{I}}$$ となり $${\mathbf{P} \neq \mathbf{S}}$$ になります。
直交回転では $${\bm \Phi = \mathbf{I}}$$ なので両者は一致します。

4.statsmodels で構造行列を得る方法

import numpy as np
import statsmodels.api as sm

# データ例
data = sm.datasets.access2.load_pandas().data.iloc[:, 1:10]

# 2因子・主軸法
fa = sm.multivariate.factor.Factor(data, n_factor=2, method='pa')
res = fa.fit()

# 斜交回転(例: promax)
res.rotate('promax')           # res.loadings がパターン行列 P になる
P = res.loadings               # (変数×因子)

# 因子間相関行列 Φ = Tᵀ T
Phi = res.rotation_matrix.T @ res.rotation_matrix

# 構造行列 S = P Φ
S = P @ Phi
  • res.loadings … パターン行列 $${\bm P}$$

  • res.rotation_matrix … 回転行列 $${\bm T}$$。斜交の場合 $${\bm \Phi = T^{\top}T}$$

  • 上記コードで S が構造行列になります。

5.まとめ

  • statsmodels が返す loadings は パターン行列。

  • 構造行列が必要なら,因子間相関を掛けて自前で算出する。

  • 直交回転なら両者は同じ,斜交回転では異なる点に注意。


2️⃣ Gemini 2.5 Flash
statsmodelsのFactor()で生成される回転後の因子負荷は、回転方法によって「因子パターン行列」と「因子構造行列」のどちらに該当するかが異なります。

  • 直交回転(例:バリマックス回転)の場合:
    直交回転では因子間に相関がないと仮定されます。
    この場合、因子パターン行列と因子構造行列は一致します。
    statsmodelsのFactorResultsオブジェクトのloadings属性で取得できるのは、この共通の因子負荷量となります。

  • 斜交回転(例:プロマックス回転、オブリミン回転)の場合:
    斜交回転では因子間に相関が許容されます。この場合、因子パターン行列と因子構造行列は異なります。

    • 因子パターン行列 (Factor Pattern Matrix): 各変数を共通因子で線形結合したときの係数(回帰係数)を表します。
      これは、項目が特定の因子から受ける純粋な影響を示します。
      解釈の際には通常、この因子パターン行列を参照します。

    • 因子構造行列 (Factor Structure Matrix): 各変数と各因子の間の相関係数を表します。
      これは、項目と因子の間の直接的な相関関係を示します。

statsmodelsのFactorResultsオブジェクトで斜交回転を行った場合、rotate()メソッドを使用すると、その結果として返されるloadingsは因子パターン行列に該当します。
因子構造行列は、因子パターン行列と因子相関行列を組み合わせて計算する必要があります。

まとめると:

  • 直交回転: loadingsは「因子パターン行列」と「因子構造行列」の両方を兼ねる。

  • 斜交回転: loadingsは「因子パターン行列」に該当する。

因子分析の結果を解釈する際には、通常、項目がどの因子に強く関連しているかを見るために因子パターン行列(斜交回転の場合)を使用することが多いです。


3️⃣ Genspark
研究結果と収集した情報に基づいて、ご質問に対する詳しい回答をまとめます。

statsmodelsの因子分析 Factor()における回転後の因子負荷について
statsmodels の因子分析(Factor())で生成される回転後の因子負荷は、因子パターン行列に該当します。

因子パターン行列と因子構造行列の違い
因子分析において、斜交回転(oblique rotation)を適用する場合、因子間に相関が許容されるため、「因子パターン行列」と「因子構造行列」という2種類の行列が生成されます。

  1. 因子パターン行列(Factor Pattern Matrix):

    • 各項目と因子の間の偏標準化回帰係数を表します

    • 他の因子の影響を制御した上での、特定の因子の各項目への固有の寄与度を示します

    • 回帰分析の偏回帰係数のように、他の因子の効果を統制した場合の影響を表します

  2. 因子構造行列(Factor Structure Matrix):

    • 各項目と因子の間の単純相関係数ゼロ次相関)を表します

    • 因子と項目の間の総合的な関連性を示します

    • 他の因子との共有効果も含めた関係を表します

直交回転と斜交回転の違い

  • 直交回転(Orthogonal Rotation):

    • 因子間の相関が0(独立)であることを前提とします

    • 直交回転の場合、因子パターン行列と因子構造行列は同一になります

    • 代表的な方法: varimax, quartimax, equamax など

  • 斜交回転(Oblique Rotation):

    • 因子間の相関を許容します

    • 因子間の相関があるため、因子パターン行列と因子構造行列は異なる値になります

    • 代表的な方法: promax, oblimin など

statsmodelsの実装について
statsmodelsのソースコードを調査した結果、以下のことが分かりました:

  1. `Factor()`クラスの`rotate()`メソッドは、回転後の因子負荷を`self.loadings`に格納します

  2. 回転には`rotate_factors()`関数が使用されています

  3. 斜交回転(promax, oblimin)でも直交回転(varimax等)でも、同様の処理で因子負荷が計算されています

  4. statsmodelsのドキュメントには、パターン行列と構造行列を明示的に区別する記述がありません

UCLAの因子分析に関するドキュメントには、「直交回転の場合、因子パターン行列と因子構造行列は同一になる」と明記されており、これは統計学の一般的な理解と一致しています。statsmodelsを含む多くの統計ソフトウェアは、回転後の因子負荷として因子パターン行列を出力することが慣例となっています。

結論
statsmodelsの`Factor()`で計算される回転後の因子負荷量(`loadings`属性)は因子パターン行列に該当します。
斜交回転を適用した場合、これは各項目と因子の間の偏回帰係数に相当し、他の因子の影響を制御した上での、特定の因子の各項目への固有の寄与度を示します。

statsmodelsには因子構造行列を直接出力する機能はないようですが、必要であれば、回転後の因子負荷(パターン行列)と因子間相関行列から別途計算することが可能です。
因子構造行列は、パターン行列に因子間相関行列を乗じることで得られます。


3つの回答は、因子パターン行列の定義と statsmodels 内部ロジックを照合していないのに、因子パターン行列と結論づけて(推測して)います

生成AIの回答を信じたいのですが…
仮説検証劇場の出力結果をどう説明したら良いのか…

本当に因子パターン行列なのかな…


シリーズの記事

次の記事

前の記事

目次

ブログの紹介


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

1.のんびり統計

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

2.実験!たのしいベイズモデリング1&2をPyMC Ver.5で

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

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

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

4.楽しい写経 ベイズ・Python等

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

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

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

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

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

7.Python機械学習プログラミング実践記

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

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

いいなと思ったら応援しよう!

ネイピア DS 応援ありがとうございます。これからもがんばって記事を作成します!

この記事が参加している募集