学生の英語テスト・数学テスト・理科テスト・国語テストの成績データがあるとします。これらの点数は互いに相関しているでしょう。しかし、「なぜ相関しているのか」を考えると、背後に「文系的な学力」と「理系的な学力」という直接観測できない潜在的な能力があり、それが各科目の成績を左右しているのではないか、と推察できます。この「直接観測できないが、観測変数の背後に存在する共通の原因」を数学的にモデル化するのが因子分析（Factor Analysis）です。
因子分析は一見すると主成分分析（PCA）と似ていますが、根本的に異なる思想に基づいています。PCAは「データの分散を最大化する方向を見つける」記述的な手法であるのに対し、因子分析は「観測データが潜在因子から生成される」という確率モデルを仮定する生成的な手法です。
因子分析を理解すると、以下のような応用が開けます。
- 心理学・教育学: 知能テストの因子構造の分析。一般知能因子（g因子）の発見は因子分析の歴史的な成果です
- マーケティング: 消費者の購買行動データから「価格志向」「品質志向」などの潜在的な因子を抽出
- 構造方程式モデリング（SEM）: 因子分析は確認的因子分析（CFA）としてSEMの測定モデルの基礎を構成します
- 遺伝学: 遺伝的変異の潜在的な要因の特定
本記事では、因子分析の数学的モデルを定式化し、パラメータの推定法を解説します。PCAとの本質的な違いを理論的に明らかにし、Pythonでの実装を通じて因子分析の動作を確認します。
本記事の内容
- 因子モデルの数学的定式化
- 共分散構造の導出と因子負荷量の意味
- PCAとの理論的な違い
- 因子回転（バリマックス回転）の理論
- 最尤推定による因子分析
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
因子モデルの数学的定式化
直感的な理解 — 「見えない原因」をモデル化する
因子分析の核心的なアイデアは、次のような問いから始まります。「観測変数の間に相関がある理由は何か？」
たとえば、身長・体重・腕の長さ・足の長さの4つの身体測定値を多くの人について測定したとします。これらの変数は互いに正の相関を持つでしょう。この相関の原因を考えると、「体のサイズ」という一つの共通因子があり、体のサイズが大きい人はすべての測定値が大きくなる、と説明できます。
因子分析はこの考え方を数学的にモデル化します。観測変数の相関を少数の共通因子で説明し、各観測変数に固有のランダムなばらつき（独自因子）を加える、という構造を仮定します。
因子モデルの定義
$p$ 個の観測変数 $\bm{x} = (x_1, x_2, \dots, x_p)^\top$ に対して、$m$ 個の共通因子 $\bm{f} = (f_1, f_2, \dots, f_m)^\top$（$m < p$）を用いた線形因子モデルは次のように定義されます。
$$ \begin{equation} \bm{x} = \bm{\mu} + \bm{\Lambda}\bm{f} + \bm{\epsilon} \end{equation} $$
ここで各記号の意味は以下の通りです。
- $\bm{x} \in \mathbb{R}^p$: 観測変数ベクトル
- $\bm{\mu} \in \mathbb{R}^p$: 観測変数の平均ベクトル
- $\bm{\Lambda} \in \mathbb{R}^{p \times m}$: 因子負荷行列（factor loading matrix）
- $\bm{f} \in \mathbb{R}^m$: 共通因子（common factors）。直接観測できない潜在変数
- $\bm{\epsilon} \in \mathbb{R}^p$: 独自因子（unique factors, specific factors）。各観測変数に固有のノイズ
因子負荷行列 $\bm{\Lambda}$ の $(i, j)$ 成分 $\lambda_{ij}$ は、第 $j$ 共通因子 $f_j$ が第 $i$ 観測変数 $x_i$ に与える影響の強さを表します。この値を因子負荷量（factor loading）と呼びます。
成分ごとに書くと、第 $i$ 観測変数は次のように表されます。
$$ x_i = \mu_i + \lambda_{i1}f_1 + \lambda_{i2}f_2 + \dots + \lambda_{im}f_m + \epsilon_i $$
これは「第 $i$ 変数の値は、$m$ 個の共通因子の線形結合に、固有のノイズ $\epsilon_i$ を加えたもの」という意味です。
モデルの仮定
因子モデルでは以下の仮定を置きます。
共通因子の仮定:
$$ \begin{equation} E[\bm{f}] = \bm{0}, \quad \text{Cov}(\bm{f}) = \bm{I}_m \end{equation} $$
共通因子の平均はゼロ、共分散行列は単位行列です。つまり、各共通因子は互いに無相関で、分散1に標準化されています。
独自因子の仮定:
$$ \begin{equation} E[\bm{\epsilon}] = \bm{0}, \quad \text{Cov}(\bm{\epsilon}) = \bm{\Psi} = \text{diag}(\psi_1, \psi_2, \dots, \psi_p) \end{equation} $$
独自因子の平均はゼロ、共分散行列 $\bm{\Psi}$ は対角行列です。つまり、異なる観測変数の独自因子は互いに無相関です。$\psi_i$ は第 $i$ 変数の独自分散（unique variance, uniqueness）と呼ばれます。
共通因子と独自因子の独立性:
$$ \begin{equation} \text{Cov}(\bm{f}, \bm{\epsilon}) = \bm{O} \end{equation} $$
共通因子と独自因子は互いに無相関（正規分布を仮定すれば独立）です。
これらの仮定は、共通因子が「純粋な潜在的原因」であり、独自因子が「各変数に固有のランダムなばらつき」であるという因子分析の思想を反映しています。
ではこれらの仮定のもとで、観測変数の共分散構造がどのように決まるか見ていきましょう。
共分散構造の導出
観測変数の共分散行列
因子モデル $\bm{x} = \bm{\mu} + \bm{\Lambda}\bm{f} + \bm{\epsilon}$ のもとで、$\bm{x}$ の共分散行列を計算します。$E[\bm{x}] = \bm{\mu}$ なので、$\bm{x} – \bm{\mu} = \bm{\Lambda}\bm{f} + \bm{\epsilon}$ です。
共分散行列の定義に従って計算します。
$$ \text{Cov}(\bm{x}) = E[(\bm{x} – \bm{\mu})(\bm{x} – \bm{\mu})^\top] $$
$\bm{x} – \bm{\mu} = \bm{\Lambda}\bm{f} + \bm{\epsilon}$ を代入すると次のようになります。
$$ \text{Cov}(\bm{x}) = E[(\bm{\Lambda}\bm{f} + \bm{\epsilon})(\bm{\Lambda}\bm{f} + \bm{\epsilon})^\top] $$
展開すると4つの項が現れます。
$$ = E[\bm{\Lambda}\bm{f}\bm{f}^\top\bm{\Lambda}^\top] + E[\bm{\Lambda}\bm{f}\bm{\epsilon}^\top] + E[\bm{\epsilon}\bm{f}^\top\bm{\Lambda}^\top] + E[\bm{\epsilon}\bm{\epsilon}^\top] $$
各項を仮定を使って整理します。$\text{Cov}(\bm{f}) = E[\bm{f}\bm{f}^\top] = \bm{I}_m$ より第1項は $\bm{\Lambda}\bm{I}_m\bm{\Lambda}^\top = \bm{\Lambda}\bm{\Lambda}^\top$ です。$\text{Cov}(\bm{f}, \bm{\epsilon}) = \bm{O}$ より第2項と第3項はゼロです。$\text{Cov}(\bm{\epsilon}) = \bm{\Psi}$ より第4項は $\bm{\Psi}$ です。
したがって、因子モデルのもとでの共分散構造は次のように表されます。
$$ \begin{equation} \bm{\Sigma} = \bm{\Lambda}\bm{\Lambda}^\top + \bm{\Psi} \end{equation} $$
これが因子分析の基本方程式です。観測変数の共分散行列 $\bm{\Sigma}$ は、因子負荷行列の積 $\bm{\Lambda}\bm{\Lambda}^\top$（共通因子による寄与）と独自分散行列 $\bm{\Psi}$（独自因子による寄与）の和に分解されます。
共通性と独自性
この分解に基づいて、第 $i$ 観測変数の分散を次のように分解できます。
$$ \sigma_{ii} = \sum_{j=1}^{m}\lambda_{ij}^2 + \psi_i $$
右辺第1項 $h_i^2 = \sum_{j=1}^{m}\lambda_{ij}^2$ を共通性（communality）と呼びます。これは、第 $i$ 変数の分散のうち共通因子で説明される部分です。
右辺第2項 $\psi_i$ は独自性（uniqueness）であり、共通因子では説明できない分散の割合です。
$$ \begin{equation} \sigma_{ii} = h_i^2 + \psi_i = \text{共通性} + \text{独自性} \end{equation} $$
変数を標準化している場合（$\sigma_{ii} = 1$）は、$h_i^2 + \psi_i = 1$ となり、共通性は「共通因子がその変数の分散の何パーセントを説明するか」を直接示します。
変数間の共分散と因子負荷
異なる変数 $i \neq j$ の共分散は次のように表されます。
$$ \sigma_{ij} = \sum_{k=1}^{m}\lambda_{ik}\lambda_{jk} $$
これは因子モデルの核心的な含意です。2つの観測変数間の相関は、共通因子を通じてのみ発生します。独自因子 $\epsilon_i$ と $\epsilon_j$ は無相関なので、共分散に寄与しません。
言い換えれば、因子分析は「変数間の相関は共通因子によって完全に説明される」と仮定しているのです。
PCAとの違いを明確にするために、次にこの2つの手法の理論的な違いを整理しましょう。
PCAとの理論的な違い
モデルの有無
PCAと因子分析の最も根本的な違いは、確率モデルの有無にあります。
PCA: データの共分散行列を固有値分解する記述的な手法です。「データが何らかの確率分布から生成された」という仮定をしません。共分散行列の固有ベクトルは常に計算でき、常に一意に定まります。PCAは「データの分散構造を要約する」だけで、データの生成メカニズムには言及しません。
因子分析: 「観測データは潜在因子から生成された」という確率モデルを仮定します。モデルのパラメータ（因子負荷行列 $\bm{\Lambda}$ と独自分散 $\bm{\Psi}$）を推定し、モデルの適合度を統計的に検定できます。因子分析は「データがなぜこのような相関構造を持つのか」を説明しようとする手法です。
共分散行列の分解の違い
この違いは、共分散行列の分解方法に如実に現れます。
PCA: $\bm{\Sigma} = \bm{W}\bm{\Lambda}_{\text{diag}}\bm{W}^\top = \sum_{k=1}^{p}\lambda_k\bm{w}_k\bm{w}_k^\top$。全分散をすべての固有値で分解します。残差項はありません。
因子分析: $\bm{\Sigma} = \bm{\Lambda}\bm{\Lambda}^\top + \bm{\Psi}$。共通因子の寄与 $\bm{\Lambda}\bm{\Lambda}^\top$ と独自分散 $\bm{\Psi}$ を明示的に分離します。対角行列 $\bm{\Psi}$ の存在が因子分析の特徴です。
独自分散 $\bm{\Psi}$ が対角行列であるという制約は、「変数間の相関は共通因子のみに起因する」という因子分析の仮定を反映しています。PCAにはこの制約がありません。
回転の不定性
因子分析には回転の不定性（rotational indeterminacy）があります。$\bm{T}$ を $m \times m$ の直交行列（回転行列）とすると、$\bm{\Lambda}$ を $\bm{\Lambda}\bm{T}$ に、$\bm{f}$ を $\bm{T}^\top\bm{f}$ に置き換えても、共分散構造は変わりません。
$$ (\bm{\Lambda}\bm{T})(\bm{\Lambda}\bm{T})^\top + \bm{\Psi} = \bm{\Lambda}\bm{T}\bm{T}^\top\bm{\Lambda}^\top + \bm{\Psi} = \bm{\Lambda}\bm{\Lambda}^\top + \bm{\Psi} = \bm{\Sigma} $$
つまり、因子負荷行列は直交回転に対して不定であり、データだけからは回転を一意に決められません。これは因子分析が「因子空間の方向」を特定できないことを意味しています。
PCAではこの不定性はありません。共分散行列の固有ベクトルは（固有値が重複しない限り）一意に決まるからです。
この回転の不定性は欠点ではなく、解釈に有利な座標系を自由に選べるという長所として活用されます。これが因子回転の考え方です。
因子回転の理論を見ていきましょう。
因子回転の理論
なぜ回転が必要か
因子分析の初期解（回転前の解）では、多くの変数が複数の因子に中程度の負荷を持つことがあり、各因子の解釈が曖昧になりがちです。理想的には、各変数は1つの因子に大きな負荷を持ち、他の因子への負荷は小さい、という単純構造（simple structure）が望まれます。
因子回転は、この単純構造に近づくように因子負荷行列を回転させる手法です。回転しても共分散構造 $\bm{\Lambda}\bm{\Lambda}^\top + \bm{\Psi}$ は変わらないので、モデルの適合度は一切損なわれません。
直交回転とバリマックス法
直交回転は直交行列 $\bm{T}$ による回転で、回転後も因子間の直交性（無相関性）が保たれます。最も広く使われる直交回転がバリマックス回転（varimax rotation）です。
バリマックス回転は、因子負荷量の2乗の分散を最大化する基準に基づいています。直感的には、「各因子に対して、負荷量が大きいか小さいかが明確に分かれる」ような状態を目指します。
具体的には、因子負荷量を標準化した上で、次のバリマックス基準 $V$ を最大化する回転行列 $\bm{T}$ を求めます。
$$ \begin{equation} V = \sum_{j=1}^{m}\left[\frac{1}{p}\sum_{i=1}^{p}\tilde{\lambda}_{ij}^4 – \left(\frac{1}{p}\sum_{i=1}^{p}\tilde{\lambda}_{ij}^2\right)^2\right] \end{equation} $$
ここで $\tilde{\lambda}_{ij} = \lambda_{ij}^*/h_i$ は共通性 $h_i$ で正規化された回転後の因子負荷量です。大括弧内は因子 $j$ における負荷量の2乗の分散に他ならず、この分散が大きいほど負荷量は「0に近い」か「大きい」かに二極化していることを意味します。
斜交回転
斜交回転（oblique rotation）は、回転後の因子が直交でなくてもよい（因子間に相関を許す）回転です。プロマックス回転（promax rotation）やオブリミン回転（oblimin rotation）が代表的です。
心理学的な因子（例えば「言語能力」と「空間認知能力」）は完全に独立ではなく、互いに相関していることが自然です。斜交回転を使えば、そのような相関のある因子構造をモデル化できます。ただし、解釈の簡便さでは直交回転に劣る場合があります。
では、因子分析のパラメータをどのように推定するか、最尤推定の理論を見ていきましょう。
パラメータの推定
主因子法
因子分析のパラメータ推定の最も古典的な方法は主因子法（principal factor method）です。この方法は反復的な手順を用います。
- 共通性の初期値を推定します（例えば、各変数の最大相関係数の絶対値の2乗を使う）
- 相関行列の対角成分を共通性で置き換えた「縮約相関行列」を計算します
- この縮約相関行列を固有値分解し、上位 $m$ 個の固有値・固有ベクトルから因子負荷行列を推定します
- 新しい因子負荷行列から共通性を再計算します
- 2〜4を共通性が収束するまで繰り返します
この方法はPCAと似ていますが、対角成分を共通性で置き換える点が本質的に異なります。PCAは全分散を扱うのに対し、主因子法は共通分散のみを扱います。
最尤推定法
因子分析の理論的に最も洗練された推定法は最尤推定法（maximum likelihood estimation, MLE）です。観測データ $\bm{x}_1, \bm{x}_2, \dots, \bm{x}_n$ が多変量正規分布 $N(\bm{\mu}, \bm{\Sigma})$ に従い、$\bm{\Sigma} = \bm{\Lambda}\bm{\Lambda}^\top + \bm{\Psi}$ という構造を持つと仮定します。
対数尤度関数は（定数項を除いて）次のように書けます。
$$ \begin{equation} \ell(\bm{\Lambda}, \bm{\Psi}) = -\frac{n}{2}\left[\ln|\bm{\Sigma}| + \text{tr}(\bm{\Sigma}^{-1}\bm{S})\right] \end{equation} $$
ここで $\bm{S}$ は標本共分散行列、$\bm{\Sigma} = \bm{\Lambda}\bm{\Lambda}^\top + \bm{\Psi}$ です。
最尤推定量は、$\ell$ を $\bm{\Lambda}$ と $\bm{\Psi}$ に関して最大化することで得られます。この最大化には解析解がなく、EM（Expectation-Maximization）アルゴリズムや勾配法などの反復アルゴリズムが使われます。
最尤推定法の大きな利点は、モデルの適合度を統計的に検定できることです。尤度比検定により、仮定した因子数 $m$ が適切かどうかを検証できます。
$$ \begin{equation} \chi^2 = n\left[\ln|\hat{\bm{\Sigma}}| – \ln|\bm{S}| + \text{tr}(\hat{\bm{\Sigma}}^{-1}\bm{S}) – p\right] \end{equation} $$
ここで $\hat{\bm{\Sigma}} = \hat{\bm{\Lambda}}\hat{\bm{\Lambda}}^\top + \hat{\bm{\Psi}}$ は推定されたモデルによる共分散行列です。この検定統計量は近似的に自由度 $\frac{1}{2}[(p-m)^2 – p – m]$ のカイ二乗分布に従います。
理論的な内容を踏まえ、Pythonで因子分析を実装してみましょう。
Pythonによる実装と可視化
合成データの生成と因子分析
まず、既知の因子構造を持つデータを生成し、因子分析がその構造を正しく復元できるかを確認します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# --- 真の因子構造 ---
n = 500 # サンプル数
p = 6 # 観測変数の数
m = 2 # 共通因子の数
# 真の因子負荷行列（2因子）
Lambda_true = np.array([
[0.9, 0.1], # x1: 因子1に強く負荷
[0.8, 0.2], # x2: 因子1に強く負荷
[0.7, 0.15], # x3: 因子1にやや強く負荷
[0.1, 0.85], # x4: 因子2に強く負荷
[0.2, 0.9], # x5: 因子2に強く負荷
[0.15, 0.75], # x6: 因子2にやや強く負荷
])
# 真の独自分散
psi_true = np.array([0.19, 0.32, 0.49, 0.26, 0.15, 0.41])
# --- データ生成 ---
F = np.random.randn(n, m) # 共通因子
E = np.random.randn(n, p) * np.sqrt(psi_true) # 独自因子
X = F @ Lambda_true.T + E
# --- 中心化と標本共分散行列 ---
X_mean = X.mean(axis=0)
X_centered = X - X_mean
S = X_centered.T @ X_centered / n
# --- 主因子法による推定 ---
# 初期の共通性推定: SMC (Squared Multiple Correlation) の近似
R = np.corrcoef(X.T) # 相関行列
np.fill_diagonal(R, 0)
h2 = np.max(np.abs(R), axis=1) # 初期共通性
for iteration in range(50):
# 縮約相関行列
R_reduced = np.corrcoef(X.T).copy()
np.fill_diagonal(R_reduced, h2)
# 固有値分解
eigvals, eigvecs = np.linalg.eigh(R_reduced)
idx = np.argsort(eigvals)[::-1]
eigvals = eigvals[idx]
eigvecs = eigvecs[:, idx]
# 上位m個の因子負荷量
Lambda_est = eigvecs[:, :m] * np.sqrt(np.maximum(eigvals[:m], 0))
# 共通性の更新
h2_new = np.sum(Lambda_est**2, axis=1)
if np.max(np.abs(h2_new - h2)) < 1e-6:
break
h2 = h2_new
# 独自分散の推定
psi_est = 1 - h2
print("=== 因子分析の結果（主因子法） ===")
print("\n因子負荷行列（推定値）:")
print(np.array2string(Lambda_est, precision=3))
print("\n因子負荷行列（真値）:")
print(np.array2string(Lambda_true, precision=3))
print(f"\n反復回数: {iteration + 1}")
この実装の結果を確認すると、推定された因子負荷行列が真の値に近い構造を持っていることがわかります。ただし、回転の不定性のため、推定された因子の方向は真の方向と一致しない場合があります。符号の反転や因子の入れ替わりが起こることもありますが、共分散構造 $\bm{\Lambda}\bm{\Lambda}^\top + \bm{\Psi}$ は同じです。
バリマックス回転の実装と可視化
次に、バリマックス回転を適用して因子の解釈可能性を高めます。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# --- データ再生成（前のコードと同じ） ---
n, p, m = 500, 6, 2
Lambda_true = np.array([
[0.9, 0.1], [0.8, 0.2], [0.7, 0.15],
[0.1, 0.85], [0.2, 0.9], [0.15, 0.75],
])
psi_true = np.array([0.19, 0.32, 0.49, 0.26, 0.15, 0.41])
F = np.random.randn(n, m)
E = np.random.randn(n, p) * np.sqrt(psi_true)
X = F @ Lambda_true.T + E
X_centered = X - X.mean(axis=0)
# 主因子法の結果（前のコードから）
R = np.corrcoef(X.T)
h2 = np.max(np.abs(R - np.eye(p)), axis=1)
for _ in range(50):
R_reduced = np.corrcoef(X.T).copy()
np.fill_diagonal(R_reduced, h2)
eigvals, eigvecs = np.linalg.eigh(R_reduced)
idx = np.argsort(eigvals)[::-1]
eigvals, eigvecs = eigvals[idx], eigvecs[:, idx]
Lambda_est = eigvecs[:, :m] * np.sqrt(np.maximum(eigvals[:m], 0))
h2_new = np.sum(Lambda_est**2, axis=1)
if np.max(np.abs(h2_new - h2)) < 1e-6:
break
h2 = h2_new
# --- バリマックス回転 ---
def varimax(Lambda, max_iter=100, tol=1e-6):
"""バリマックス回転の実装"""
p, m = Lambda.shape
R_mat = np.eye(m)
L = Lambda.copy()
for _ in range(max_iter):
old_L = L.copy()
for i in range(m):
for j in range(i+1, m):
# 2因子ずつ回転角を最適化
u = L[:, i]**2 - L[:, j]**2
v = 2 * L[:, i] * L[:, j]
A = np.sum(u)
B = np.sum(v)
C = np.sum(u**2 - v**2)
D = 2 * np.sum(u * v)
# 回転角
num = D - 2 * A * B / p
den = C - (A**2 - B**2) / p
phi = 0.25 * np.arctan2(num, den)
# 回転の適用
cos_phi = np.cos(phi)
sin_phi = np.sin(phi)
L_new_i = cos_phi * L[:, i] + sin_phi * L[:, j]
L_new_j = -sin_phi * L[:, i] + cos_phi * L[:, j]
L[:, i] = L_new_i
L[:, j] = L_new_j
if np.max(np.abs(L - old_L)) < tol:
break
return L
Lambda_rotated = varimax(Lambda_est)
# --- 可視化 ---
var_names = [f"$x_{i+1}$" for i in range(p)]
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) 回転前の因子負荷量
ax = axes[0]
im = ax.imshow(Lambda_est, cmap="RdBu_r", vmin=-1, vmax=1, aspect="auto")
ax.set_xticks(range(m))
ax.set_xticklabels(["Factor 1", "Factor 2"], fontsize=10)
ax.set_yticks(range(p))
ax.set_yticklabels(var_names, fontsize=11)
ax.set_title("(a) Before rotation", fontsize=12)
for i in range(p):
for j in range(m):
ax.text(j, i, f"{Lambda_est[i,j]:.2f}", ha="center", va="center",
fontsize=10, color="white" if abs(Lambda_est[i,j]) > 0.5 else "black")
plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
# (b) 回転後の因子負荷量
ax = axes[1]
im = ax.imshow(Lambda_rotated, cmap="RdBu_r", vmin=-1, vmax=1, aspect="auto")
ax.set_xticks(range(m))
ax.set_xticklabels(["Factor 1", "Factor 2"], fontsize=10)
ax.set_yticks(range(p))
ax.set_yticklabels(var_names, fontsize=11)
ax.set_title("(b) After varimax rotation", fontsize=12)
for i in range(p):
for j in range(m):
ax.text(j, i, f"{Lambda_rotated[i,j]:.2f}", ha="center", va="center",
fontsize=10, color="white" if abs(Lambda_rotated[i,j]) > 0.5 else "black")
plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
# (c) 因子負荷量の散布図
ax = axes[2]
ax.scatter(Lambda_est[:, 0], Lambda_est[:, 1], c="blue", s=80,
alpha=0.7, label="Before rotation", marker="o")
ax.scatter(Lambda_rotated[:, 0], Lambda_rotated[:, 1], c="red", s=80,
alpha=0.7, label="After rotation", marker="^")
for i in range(p):
ax.annotate(f"$x_{i+1}$",
(Lambda_rotated[i, 0] + 0.03, Lambda_rotated[i, 1] + 0.03),
fontsize=10, color="red")
ax.set_xlabel("Factor 1 loading", fontsize=11)
ax.set_ylabel("Factor 2 loading", fontsize=11)
ax.set_title("(c) Loading scatter plot", fontsize=12)
ax.axhline(0, color="k", linewidth=0.5)
ax.axvline(0, color="k", linewidth=0.5)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(-0.2, 1.1)
ax.set_ylim(-0.2, 1.1)
ax.set_aspect("equal")
plt.tight_layout()
plt.savefig("factor_analysis_rotation.png", dpi=150, bbox_inches="tight")
plt.show()
この可視化から、因子回転の効果が明確に読み取れます。
-
回転前（図a）: 多くの変数が両方の因子に中程度の負荷を持っており、「因子1は何を表すのか」「因子2は何を表すのか」の解釈が曖昧です
-
回転後（図b）: バリマックス回転を適用すると、$x_1, x_2, x_3$ は因子1に強く負荷し因子2にはほとんど負荷せず、逆に $x_4, x_5, x_6$ は因子2に強く負荷し因子1にはほとんど負荷しません。真の因子構造に近い「単純構造」が得られています
-
因子負荷量の散布図（図c）: 回転前は変数が2つの因子の中間方向に散らばっていますが、回転後は軸に近い位置に集まっています。これが単純構造の幾何学的な意味です
scikit-learnを用いた因子分析
実用的な因子分析にはscikit-learnの FactorAnalysis クラスが便利です。最尤推定法に基づくパラメータ推定が実装されています。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import FactorAnalysis, PCA
from sklearn.datasets import load_iris
# --- データの読み込み ---
iris = load_iris()
X = iris.data
y = iris.target
feature_names = ["sepal_len", "sepal_wid", "petal_len", "petal_wid"]
# --- 因子分析（2因子） ---
fa = FactorAnalysis(n_components=2, random_state=42)
F_scores = fa.fit_transform(X)
Lambda_fa = fa.components_.T # (p, m)
psi_fa = fa.noise_variance_
# --- PCA（2成分） ---
pca = PCA(n_components=2)
PC_scores = pca.fit_transform(X)
Lambda_pca = pca.components_.T * np.sqrt(pca.explained_variance_)
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) PCA vs FA のスコアプロット
colors = ["#e41a1c", "#377eb8", "#4daf4a"]
ax = axes[0]
for k in range(3):
mask = y == k
ax.scatter(F_scores[mask, 0], F_scores[mask, 1],
c=colors[k], alpha=0.6, s=25, label=iris.target_names[k])
ax.set_xlabel("Factor 1", fontsize=11)
ax.set_ylabel("Factor 2", fontsize=11)
ax.set_title("(a) Factor Analysis scores", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (b) 因子負荷量の比較
ax = axes[1]
x_pos = np.arange(len(feature_names))
width = 0.35
ax.barh(x_pos - width/2, Lambda_fa[:, 0], width, color="steelblue",
alpha=0.7, label="FA Factor 1")
ax.barh(x_pos + width/2, Lambda_fa[:, 1], width, color="coral",
alpha=0.7, label="FA Factor 2")
ax.set_yticks(x_pos)
ax.set_yticklabels(feature_names, fontsize=10)
ax.set_xlabel("Factor loading", fontsize=11)
ax.set_title("(b) Factor loadings (FA)", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, axis="x")
ax.axvline(0, color="k", linewidth=0.5)
# (c) 独自分散
ax = axes[2]
ax.barh(x_pos, psi_fa, color="gray", alpha=0.7)
ax.set_yticks(x_pos)
ax.set_yticklabels(feature_names, fontsize=10)
ax.set_xlabel("Unique variance $\\psi_i$", fontsize=11)
ax.set_title("(c) Uniqueness", fontsize=12)
ax.grid(True, alpha=0.3, axis="x")
plt.tight_layout()
plt.savefig("factor_analysis_iris.png", dpi=150, bbox_inches="tight")
plt.show()
print("=== 因子負荷行列 ===")
for i, name in enumerate(feature_names):
print(f"{name}: F1={Lambda_fa[i,0]:.3f}, F2={Lambda_fa[i,1]:.3f}, "
f"h²={np.sum(Lambda_fa[i]**2):.3f}, ψ={psi_fa[i]:.3f}")
print(f"\n共分散再構成の相対誤差: "
f"{np.linalg.norm(fa.get_covariance() - np.cov(X.T)) / np.linalg.norm(np.cov(X.T)):.4f}")
アイリスデータに対する因子分析の結果から、以下の知見が得られます。
-
因子スコアプロット（図a）: PCAのスコアプロットと似た構造が得られますが、微妙に異なっています。因子分析は独自分散を明示的にモデル化しているため、独自分散が大きい変数の影響が相対的に抑えられ、共通因子による変動がより強調されます
-
因子負荷量（図b）: petal_lenとpetal_widが因子1に大きく負荷しており、sepal_widが因子2に負荷しています。これはPCAの因子負荷量と同じ傾向ですが、因子分析では独自分散を分離しているため、負荷量の値が異なります
-
独自分散（図c）: 各変数の独自分散が推定されています。独自分散が大きい変数は、共通因子ではうまく説明できない固有の変動を持っていることを意味します。sepal_widの独自分散が比較的大きいのは、この変数が他の変数とあまり相関しないことの反映です
まとめ
本記事では、因子分析の理論を数学的に定式化し、その特徴を解説しました。
- 因子モデル: $\bm{x} = \bm{\mu} + \bm{\Lambda}\bm{f} + \bm{\epsilon}$ は、観測変数を少数の共通因子と独自因子に分解する確率モデルです。共分散構造は $\bm{\Sigma} = \bm{\Lambda}\bm{\Lambda}^\top + \bm{\Psi}$ と表されます
- PCAとの違い: PCAはデータの分散構造を記述する手法で、因子分析はデータの生成メカニズムを仮定するモデルベースの手法です。因子分析は独自分散 $\bm{\Psi}$ を明示的にモデル化する点が本質的に異なります
- 因子回転: 共分散構造を変えずに因子負荷行列を回転させることで、解釈しやすい単純構造を得ることができます。バリマックス回転は最も広く使われる直交回転法です
- パラメータ推定: 主因子法と最尤推定法が代表的です。最尤推定法はモデルの適合度検定が可能であり、理論的に好ましい推定法です
因子分析は、データの背後にある潜在的な構造を探索する強力な道具ですが、因子数の決定や回転方法の選択など、分析者の判断に委ねられる部分も多い手法です。結果の解釈には専門知識と慎重な検討が必要です。
次のステップとして、以下の記事も参考にしてください。
- 主成分分析の理論 — PCAの理論をより深く理解する
- 正準相関分析 — 2組の変数群間の関係を分析する手法
- パス解析と構造方程式モデリング — 因子分析を発展させた因果モデリング