Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
8
Help us understand the problem. What are the problem?

posted at

updated at

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

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

因子分析の概要

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

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

fa.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です。

変数名 意味
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の方は回転なども扱えるようだったので、今回はこちらを使用します。

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

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に格納します。

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]

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

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 を指定します。

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

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

In
fa.loadings_
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() を呼び出します。可視化のためにラベルも加えてみます。

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

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

各因子の考察

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

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

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

  • 【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_phenolsflavanoids といった値に影響を与えているのであれば、分析には意味があると言えます。化学に知見があればもう少し深い考察ができるかもしれません。

因子スコアによる分類

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

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)
In
acc = accuracy_score(y, y_pred)
# 0.9831460674157303

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

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

参考

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
8
Help us understand the problem. What are the problem?