Python scikit-learnによる因子分析

     本記事ではsklearnのFactorAnalysisを利用した。因子分析は主成分分析と同様に説明変数Xとは別の変数を定義して分析する手法である。結果的に次元削減を目的とする。

　因子分析は既存の説明変数Xを分解して、「共通因子」なる仮定(想定)の変数での構成する次式(1)に変換する。共通因子Factor1,2,…に対してそれぞれの負荷量Loadとの積とし、最後に独自因子と呼ぶ誤差eiを加えている。
Xi=Factor1×Load1＋Factor2×Load2+ … +ei …(1)
　特記事項は、共通因子の数は自分で決めることである。

  データ分析には、sklearnに含まれるデータセットbostonを使用した。目的変数Yは住宅価格で、説明変数Xは部屋数や犯罪率など13因子ある。因子分析した結果を下図4つを示す。想定する共通因子を2/3/12/13因子としてそれぞれ分析した場合である。各図の横軸は共通因子Factor1、縦軸は共通因子Factor2として全データをプロットしている。

　4枚の図は見た目に差がある。各図の横軸はFactor1で縦軸はFactor2であるが4枚の図間で同じ意味でない、これは式(1)よりわかる。

    各図の見方について説明する。図中の赤矢印の向きと大きさは共通因子に対する各説明変数Xiの影響度を表す。例えば、方向が同じ説明変数Xiが複数あった場合、互いに依存性をもっている可能性もあり得る。他には、赤矢印の方向に生データプロットが少なければ、それらを外れ値成分と判断することも考える。これにより、幾つかの説明変数Xiを省くことによる次元削減の検討を出来はする。

▼本プログラム

# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from sklearn import datasets
from sklearn import linear_model
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import FactorAnalysis as FA

def my_load_dataset(): # データセットを読み込みXとYで返す
    boston = datasets.load_boston()
    boston_df=pd.DataFrame(boston.data) # 説明変数Xをpandas形式に変換
    boston_df.columns = boston.feature_names # カラム名を挿入
    boston_df['PRICE'] = pd.DataFrame(boston.target) # 目的変数Yをカラム名PRICEで追記   
    x_df = boston_df.drop("PRICE", axis=1) # Xデータのみ抽出。pandas形式で列名有り
    #Y_df = boston_df.loc[:,['PRICE']] # Yデータのみ抽出。pandas形式で列名有り
    y = boston_df.PRICE # Yデータのみ抽出。カラム名を指定して抽出のため、列名はない
    df_columns = boston.feature_names
    return x_df, y, df_columns # XYデータをタプルで返す。[]を付けた場合はリストとなる。

def my_scaler(X, df_columns): # 説明変数Xを標準化
    scaler = StandardScaler()
    x_sc = scaler.fit_transform(X)
    ### plot    
    # sklearnの上記のようなfit関連の実行後、データの書式がnumpy形式となりカラム(列)名が無くなる。
    # そのため、再度、pandas形式に変換する(次の2行)。これはグラフの軸名で使用するため。
    x_sc = pd.DataFrame(x_sc)
    x_sc.columns = df_columns
    return x_sc

def my_fa(X_sc, Y, n): # 因子分析
    fa = FA(n_components=n, max_iter=500)
    fa.fit(X_sc)
    print("fa.components_")
    print(fa.components_)
    
    x_fa = fa.fit_transform(X_sc)
    ### plot
    x_data = x_fa[:,0]
    y_data = x_fa[:,1]
    fa0 = fa.components_[0]
    fa1 = fa.components_[1]
    plt.figure(figsize=(9,9))
    plt.scatter(x_data, y_data, c=Y/len(set(Y)), marker=".")
    for i in range(fa0.shape[0]):
        plt.arrow(0, 0, fa0[i]*1.5, fa1[i]*1.5, color='r', width=0.0005, head_width=0.07, alpha=0.8)
        plt.text(fa0[i]*1.8, fa1[i]*1.8, X_sc.columns.values[i], color='r', fontsize = 18)
    plt.xlim([-2.5,2.5])
    plt.ylim([-2.5,2.5])
    plt.xlabel('Factor 1', fontsize = 20)
    plt.ylabel('Factor 2', fontsize = 20)
    plt.title(str(n)+" factor", fontsize = 20)
    plt.tick_params(labelsize = 18)
    plt.grid()
    c = patches.Circle(xy=(0, 0), radius=1, ec='r', fill=False)
    plt.axes().add_patch(c)
    #plt.show()
    plt.savefig("FA_" + str('%02g' %n) + '.png')

def my_regression(X, Y, mytitle='Ref', mylabel='Multiple regression analysis'): # 線形重回帰分析
    clf = linear_model.LinearRegression(fit_intercept=True, normalize=False, copy_X=True, n_jobs=1)
    clf.fit(X, Y)
    print(pd.DataFrame({"Name":X.columns, "Coefficients":clf.coef_}).sort_values(by='Coefficients')) # 係数でソート
    print(clf.intercept_) # 切片
    print(clf.score(X,Y)) # R^2

    plt.figure()
    plt.title(mytitle)
    plt.xlabel('RM (number of rooms)', fontsize=14)
    plt.ylabel('PRICE (target)', fontsize=14)
    plt.scatter(X["RM"], Y, c='blue', label='Raw data')
    plt.scatter(X.RM, clf.predict(X), c='red', label=mylabel)
    plt.legend(loc='lower right', fontsize=12)
    plt.show()

if __name__ == '__main__':
    # データセットを読み込み次のの3つを返す。X, Y, 説明変数Xのカラム(列)名
    data = my_load_dataset()
    X = data[0]
    Y = data[1]
    df_columns = data[2]
    print(X)
    
    # 説明変数Xを標準化
    X_sc = my_scaler(X, df_columns)
    print(X_sc)
    
    # 因子分析
    mycheck_Num=[2,3,12,13]
    for i in mycheck_Num:
        print(i)
        my_fa(X_sc, Y, i)
    
    # 説明変数Xを減らす
    # 　因子分析の結果から独自判断で下記4つの説明変数を削除
    drop_col_1 = ['LSTAT', 'INDUS', 'NOX', 'AGE']
    X1 = X.drop(drop_col_1, axis=1) # X1は9つの説明変数Xで構成されるデータフレーム(pandas形式)

    # 線形重回帰分析
    my_regression(X, Y, mytitle='Ref', mylabel='13feature analysis')  # 1.説明変数X全13因子を使用しての回帰
    my_regression(X1, Y, mytitle='After FA', mylabel='9feature analysis') # 2.説明変数X9因子を使用しての回帰

参考リンク

sklearn.decomposition.FactorAnalysis — scikit-learn 0.19.1 documentation

Python scikit-learnによる主成分分析 - HK29’s blog

    参考までに、主成分分析と因子分析の違いを記載する。主成分分析の場合は、目的変数Yに対する説明変数を主成分PCi＝X1×e1＋X2×e2＋…と呼ぶ合成変数に置換している。その定義は分散に着目しており、第一主成分PC1が最大の分散方向に軸をとった成分で、最も寄与率が高いため、分析時のコンポーネント数が2だろうが、10だろうがPC1の結果に違いはない。
　一方、因子分析の場合は、式(1)をみるとわかる通り共通因子Factorは、説明変数Xiに従属しており、目的変数Yに対する説明変数はXiのままである。そのため、共通因子Factorの数に依って負荷量Loadの感度(大きさ)は異なる。別の表現をすると、その数に依存して相対関係に変化は生じず、同じ解釈に成りうる。結局は、人の判断に依る。

    データ分析手法は色々あるが、最終目的は回帰式の作成と利用であろうから、個人的にはいきなり重回帰分析、もしくはランダムフォレストで良いと思う。

以上

<広告>

リンク