Pythonプログラミング（ステップ７・統計計算・因子分析）

このページでは、NumPyとRを使った因子分析について考える。

対数尤度の最大化

サンプルの共分散行列$C$がわかったときの対数尤度 $$ L(C; W, \Psi) = -\frac{n N}{2} \log(2 \pi) - \frac{n}{2} \log\left| WW^\top + \Psi \right| - \frac{n}{2} \textrm{tr} \left( (W W^\top + \Psi)^{-1} C \right) $$ を最大化する方法を考えよう。

天下りではあるが、データとモデルの分散・共分散が一致する場合、すなわち、 $$ W W^\top + \Psi = C \tag{1} $$ のとき、対数尤度は、$C$の固有値$\lambda_i$を使って $$ L(C; W, \Psi) = -\frac{n N}{2} \log(2 \pi) - \frac{n}{2} \sum_i^N \log \lambda_i - \frac{nN}{2} $$ となり、これが尤度の最大値を与える。

ここで、問題なのは、このままでは$W$と$\Psi$の中の未知数のほうが$C$の独立な成分の数よりも多くなると、これだけの情報からは未知数を決めることができない（不良設定問題になる）ことである。また、仮定する独立な「ノイズ源」の数（共通因子）が少ないと、反対に、(1)式を完全に満たすような$W$と$\Psi$は存在しなくなってしまう。

対数尤度が最大となる条件

統計パッケージの活用

以上のように、モデルのパラメータ（$W$と$\Psi$）、および因子の数($W$の列の数）を適切に決めるのは、 様々な配慮すべき事項があって、それを処理するコードを作成するにも職人技が必要となる。 そこで、専門家によって統計処理用の専用のソフトウェアが開発され、広く普及している。 ここでは、代表的な統計ソフトであるRの機能をPythonから呼び出し、因子分析を行わせてみよう。

PythonからRの機能を利用する手段はいくつか提供されているが、ここでは rpy2を用いることにする。 rpy2はPythonの標準機能ではないので、追加インストールする必要がある。 また、rpy2は内部でRを呼び出すので、Rも動作する状態にセットアップしておく必要がある。

以下は、主成分分析のページで用いた、プロ野球球団選手のサンプルデータ （身長、体重、年齢、在籍年数、年収）に対して因子分析を行うコードの例である：

# coding: utf-8

import matplotlib.pyplot as plt
from rpy2.robjects import r, numpy2ri
import numpy as np

numpy2ri.activate() 

data = np.loadtxt("rakuten-2018.txt", delimiter=" ", skiprows=1, usecols=(0,1,2,3,4) )
ndata = data.shape[0]

data[:,4] = np.log10(data[:,4])

factanal = r('factanal')
res = factanal(data,2,'promax')

print(res[1])

print('Uniqueness:')
print(res[2])

corr_mat = res[3]
eig = np.linalg.eig(corr_mat)
print('Eigenvalues:')
print(eig[0])

X = np.array(res[1])
plt.scatter(X[:,0],X[:,1])
labels=['height','weight','age','years','salary']
for i,label in enumerate(labels):
    plt.annotate(label, (X[i,0], X[i,1]))
plt.grid()
plt.xlabel('Factor 1')
plt.xlim(-1,1)
plt.ylabel('Factor 2')
plt.ylim(-1,1)
plt.show()

コードの説明

from rpy2.robjects import r, numpy2ri

rpy2の機能のうち、Rオブジェクトの操作を r(何々)またはr.何々、RとNumPyオブジェクトとの変換機能を numpy2ri.何々 で利用できるようにする。

numpy2ri.activate()

RとNumPyオブジェクトとの変換機能を「活性化」させる。

factanal = r('factanal')

Rの factanal（因子分析を行うための機能）をPythonコードの factanal に対応づける。

res = factanal(data,2,"promax")

NumPyのオブジェクトdataに対して、共通因子の数を 2 とし、因子の回転に promax 法を用い、因子分析を行う（その他、"none", "varimax" が選択できる）。 ここで、データは自動的に標準化（平均0、分散1に調整）された上で処理される。 結果を res に代入。その際、Rの出力オブジェクトの要素は、Python(rpy2)側ではリストに格納される。

print(res[1])

因子分析の結果のうち、各変数の因子への負荷を出力する。

print(res[2])

因子分析の結果のうち、各項目の独自性を出力する。

Scikit-learnを用いる例

# coding: utf-8

import matplotlib.pyplot as plt
from sklearn.decomposition import FactorAnalysis
import numpy as np

data = np.loadtxt("rakuten-2018.txt", delimiter=" ", skiprows=1, usecols=(0,1,2,3,4) )
ndata = data.shape[0]

data[:,4] = np.log10(data[:,4])

mean = np.mean(data,0)
stdev = np.std(data,0)

for i in range(5):
    data[:,i] = (data[:,i]-mean[i])/stdev[i]

fa = FactorAnalysis(2)
fa.fit(data)

corr_mat = fa.get_covariance()
eig = np.linalg.eig(corr_mat)
print('Eigenvalues:')
print(eig[0])

X = fa.components_.T
print('Loadings:')
print(X)

plt.scatter(X[:,0],X[:,1])
labels=['height','weight','age','years','salary']
for i,label in enumerate(labels):
    plt.annotate(label, (X[i,0], X[i,1]))
plt.grid()
plt.xlabel('Factor 1')
plt.xlim(-1,1)
plt.ylabel('Factor 2')
plt.ylim(-1,1)
plt.show()

　練習：球団の因子分析

上のPythonコードの動作を確認しなさい。Rのfactanalの動作について調べた上で、共通因子の数、回転の方法（有無）を変えつつ、 出力結果を解釈しなさい。

　解説:「回転」操作

因子分析で前提としているモデル $$ X = W F + \Psi $$ において、座標変換を表す行列 $Q$ を使い、 $$ X = W Q Q^{-1} F + \Psi $$ とおいて、 $$ W \leftarrow W Q $$ かつ $$ F \leftarrow Q^{-1} F $$ であったと思い直しても、話の本質は何も変わらない（変換して逆戻しするだけだから）。 もし、$Q$が回転を表す行列（直交行列）であれば$Q^{-1} = Q^\top$である。 そこで、結果の解釈が分かりやすい（因子負荷量が「綺麗に」グループ化される）ように、うまい具合に$Q$を選んで「前処理」しておくと都合が良い。 factanal( )の「因子の回転」のオプションは、こうした動作を指定するために設けられている。

　練習：旅行サイトのレビューの因子分析

旅行サイトのレビューの実データ（CSV形式）を こちら (東北大学限定)に用意した。

データは、「施設ID」, 「立地」, 「部屋」, 「食事」, 「風呂」, 「サービス」, 「設備」, 「総合」 の項目で構成され、 「施設ID」以外の、「立地」から「総合」までの７項目についてそれぞれ1点から5点までの5段階で点数化されている。

上のサンプルコードを手直しの上、このデータを因子分析し、結果を解釈してみなさい。

　ヒント

data = np.loadtxt('travel_data.csv',delimiter=',',dtype='float32',skiprows=1,usecols=(1,2,3,4,5,6,7))