Pythonで因子分析 〜特徴量から意味のある因子を抽出する〜

因子分析を概説しつつ、Pythonで動かしてみよう記事です。

## 因子分析の概要

多数の特徴量からなるデータから、特徴量間の共通の因子を探り出し、少数の共通因子によってデータを単純化する手法を**因子分析** (factor analysis) といいます。

5教科のテストをn人の生徒が受験した時、各科目の点数はそれぞれ独立に決定されるわけではなく、例えば「理系の能力」が高い人は「数学」と「理科」の点数がともに高い傾向にあり、「文系の能力」が高い人は「国語」「英語」「社会」が高くなる傾向にあると考えられます。このように「理系の能力」「文系の能力」といった共通の要因を**共通因子** (common factor) といいます。

![fa.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/fa3a4f23-09ef-b9bc-0ee5-4c4560a25eed.png)


「理系の能力」「文系の能力」といった因子の各生徒が持つ値を**因子スコア** (factor score) といい、今、生徒iのk番目の因子スコアを$f_{ik}$と表すことにします。このとき、K個の因子の存在を仮定するモデル (K因子モデル) において、(生徒i, 科目j) の点数$x_{ij}$は以下のように表されます。

$$x_{ij}=a_{j1}f_{i1}+...+a_{jK}f_{iK}+d_{j}u{ij}$$

ここで$a_{jk}$は**因子負荷量**と呼ばれ、因子分析で推定される対象になります。また$d_{j}u{ij}$の項は共通因子では説明できなかった部分を表しており、$u{ij}$は独自性因子、$d_{j}$は独自係数などと呼ばれます。
上の式を行列で表現すると以下のようになります。

$$X=FA^{\mathrm{T}}+UD$$

因子負荷量を要素に持つ$A$を**因子負荷行列**といいます。
なお、因子の解釈が分かりやすくなるような座標軸を探索する**回転** (rotation) という操作がありますが、本記事では触れません。有名なものとしては**バリマックス回転** (varimax rotation) などがあります。

それでは実際にデータを使って因子分析を試してみましょう！

# 使用するデータセットについて

イタリアの同じ地区で生産された3つの異なる品種のワインの化学分析結果に関するデータです。説明変数として以下の13個の特徴量があり、目的変数としてワインの品種 (`class_0`, `class_1`, `class_2`) が用意されています。データサイズは178です。

* [Wine Data Set - UCI Machine Learning Repository](https://archive.ics.uci.edu/ml/datasets/wine)

|変数名|意味|
|:--|:--|
|alcohol|アルコール度数|
|malic_acid|リンゴ酸|
|ash|灰 (？)|
|alcalinity_of_ash|灰 (？) のアルカリ度|
|magnesium|マグネシウム|
|total_phenols|総フェノール類量|
|flavanoids|フラボノイド (ポリフェノールの一種)|
|nonflavanoid_phenols|非フラボノイドフェノール類|
|proanthocyanins|プロアントシアニジン (ポリフェノールの一種)|
|color_intensity|色の強さ|
|hue|色合い|
|od280/od315_of_diluted_wines|ワインの希釈度合い|
|proline|プロリン (アミノ酸の一種)|

# factor_analyzerによる因子分析

因子分析のためのPythonパッケージとしては主に次の2つがあるようです。

* `sklearn.decomposition.FactorAnalysis`
* `factor_analyzer`

factor_analyzerの方は回転なども扱えるようだったので、今回はこちらを使用します。

## 必要なライブラリのインポート

```python:In
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.datasets import load_wine
from factor_analyzer import FactorAnalyzer
```

## データの読込み・前処理

scikit-learnが提供しているデータセットを読み込み、pandas DataFrameに格納します。

```python:In
wine = load_wine()
X = pd.DataFrame(wine.data, columns=wine.feature_names)
print(X)
```

```:Out

     alcohol  malic_acid   ash  alcalinity_of_ash  magnesium  total_phenols  flavanoids  nonflavanoid_phenols  proanthocyanins  color_intensity   hue  od280/od315_of_diluted_wines  proline
0      14.23        1.71  2.43               15.6      127.0           2.80        3.06                  0.28             2.29             5.64  1.04                          3.92   1065.0
1      13.20        1.78  2.14               11.2      100.0           2.65        2.76                  0.26             1.28             4.38  1.05                          3.40   1050.0
2      13.16        2.36  2.67               18.6      101.0           2.80        3.24                  0.30             2.81             5.68  1.03                          3.17   1185.0
3      14.37        1.95  2.50               16.8      113.0           3.85        3.49                  0.24             2.18             7.80  0.86                          3.45   1480.0
4      13.24        2.59  2.87               21.0      118.0           2.80        2.69                  0.39             1.82             4.32  1.04                          2.93    735.0
..       ...         ...   ...                ...        ...            ...         ...                   ...              ...              ...   ...                           ...      ...
173    13.71        5.65  2.45               20.5       95.0           1.68        0.61                  0.52             1.06             7.70  0.64                          1.74    740.0
174    13.40        3.91  2.48               23.0      102.0           1.80        0.75                  0.43             1.41             7.30  0.70                          1.56    750.0
175    13.27        4.28  2.26               20.0      120.0           1.59        0.69                  0.43             1.35            10.20  0.59                          1.56    835.0
176    13.17        2.59  2.37               20.0      120.0           1.65        0.68                  0.53             1.46             9.30  0.60                          1.62    840.0
177    14.13        4.10  2.74               24.5       96.0           2.05        0.76                  0.56             1.35             9.20  0.61                          1.60    560.0

[178 rows x 13 columns]
```

値の大小にばらつきがあるので正規化を行います。

```python:In
X_stats = X.describe().transpose()

def norm(x):
    return (x-X_stats['mean'])/X_stats['std']

X = norm(X)
```

## 因子分析

モデルを生成し、`fit()` メソッドでデータに適合させます。`n_factors` には仮定する因子の数 (default:3) を設定し、回転を与える場合には `rotation` を指定します。

```python:In
fa = FactorAnalyzer(n_factors=3, rotation=None)
fa.fit(X)
```

計算された因子負荷行列は `loadings_` 属性に格納されています。

```python:In
fa.loadings_
```

```r:Out
array([[ 3.03695933e-01,  7.18459286e-01, -1.84135256e-01],
       [-4.69444771e-01,  2.78685460e-01,  8.61921144e-02],
       [-8.96586919e-04,  4.49071087e-01,  6.51789626e-01],
       [-5.10072131e-01, -2.83591087e-02,  6.99150855e-01],
       [ 2.63599880e-01,  3.49076632e-01,  9.10060413e-02],
       [ 8.45207510e-01,  9.47272285e-02,  1.66562978e-01],
       [ 9.38545127e-01, -1.21792372e-02,  1.90407577e-01],
       [-5.80889756e-01,  4.43611005e-02,  1.38685529e-01],
       [ 6.18386558e-01,  4.26912975e-02,  1.22035991e-01],
       [-1.85667896e-01,  8.09678094e-01, -1.08117369e-01],
       [ 5.98261734e-01, -3.88603476e-01,  5.02801023e-02],
       [ 8.03486817e-01, -2.62017974e-01,  1.74507345e-01],
       [ 6.02834445e-01,  5.48470217e-01, -1.24683648e-01]])
```

観測データを因子負荷行列を使って因子スコアに変換するには、引数にデータを与えて `transform()` を呼び出します。可視化のためにラベルも加えてみます。

```python:In
result_df = pd.DataFrame(fa.transform(X), columns=['factor1', 'factor2', 'factor3'])
result_df['class'] = wine.target

sns.scatterplot(x='factor1', y='factor2', hue='class', data=result_df)
```

![スクリーンショット 2020-12-05 12.10.19.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/4ba5b131-f5f2-004a-a54e-b81d1791751c.png)

factor1とfactor2を軸にとってみたところ、何となくクラス毎に別れている様子が確認できます。

## 各因子の考察

さて、因子分析によって導かれた各因子はどのような意味を持つでしょうか。見やすいように因子負荷行列をDataFrameに入れてみます。

```python:In
loadings_df = pd.DataFrame(fa.loadings_, columns=['factor1', 'factor2', 'factor3'])
loadings_df.index = wine.feature_names

loadings_df
```

![スクリーンショット 2020-12-05 12.13.51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/a93d70df-c647-f722-1426-03b8f70752ad.png)


特徴量が専門的なので解釈が難しいですが、以下の傾向が見えるような気がします。

* **【factor1】** :  
    - (>0.6) `total_phenols`, `flavanoids`, `proanthocyanins`, `od280/od315_of_diluted_wines`, `proline	`
    - フェノール関係の特徴量が集まってるっぽい？
* **【factor2】** : 
    - (>0.6) `alcohol`, `color_intensity`
    - アルコール度数と色の強さ (関連性はよく分からない) 
* **【factor3】** : 
    - (>0.6) `ash`, `alcalinity_of_ash`
    - 灰 (?) に関する特徴量

何らか共通の因子 (例えばfactor1) が存在して、それが原因で `total_phenols` や `flavanoids` といった値に影響を与えているのであれば、分析には意味があると言えます。化学に知見があればもう少し深い考察ができるかもしれません。

# 因子スコアによる分類

せっかくなので分類までやってみましょう。得られた因子スコアのデータとラベルを使ってSVMで分類してみます。

```python:In
from sklearn.model_selection import KFold
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score, confusion_matrix

svc = SVC(C=1.0)
y = wine.target

# 5分割交差検証
kf = KFold(n_splits=5, random_state=0, shuffle=True)

# 予測結果格納用NumPy配列
y_pred = np.zeros((len(X),), dtype=float)

for train_idx, valid_idx in kf.split(X):
    X_train, X_valid = X.iloc[train_idx], X.iloc[valid_idx]
    y_train, y_valid = y[train_idx], y[valid_idx]
    svc.fit(X_train, y_train)
    y_pred[valid_idx] = svc.predict(X_valid)
```

```python:In
acc = accuracy_score(y, y_pred)
# 0.9831460674157303

cnf = confusion_matrix(y, y_pred)
# [[58  1  0]
#  [ 0 70  1]
#  [ 0  1 47]]
```

次元数を落とした上でしっかり分類できています。めでたし！

# 参考

* [日本統計学会公式認定 統計検定準1級対応 統計学実践ワークブック](https://www.amazon.co.jp/dp/478060852X/)
* [factor_analyzer Documentation Release 0.3.1](https://buildmedia.readthedocs.org/media/pdf/factor-analyzer/latest/factor-analyzer.pdf)
* [2. Pythonで綴る多変量解析　4-1. 因子分析(scikit-learn) - Qiita](https://qiita.com/y_itoh/items/227cb33317ceb09199c2)