Python scikit-learnによる主成分分析

　本記事では、sklearnのPrincipal component analysis (PCA)を利用した。主成分分析の目的は次元削減であり、説明変数Xそのものではなく分散に着目した分析方法である。はじめに、全データの分散が最大となる方向に軸をとって第一主成分と呼ぶ。これは複数の説明変数Xから合成変数として新規に定義したことを意味する。その第一主成分に与える影響度の低い説明変数Xを調べることで次元削減を目的とする。更に、第一主成分に垂直な方向の軸を第二主成分と呼び、同様の操作を順に行う。

　主成分分析での注意点は、説明変数Xを標準化する必要がある。全分散に着目した分析方法のため、各説明変数Xの平均を0にして標準偏差を1にして説明変数間を平準化する。これはsklearnのStandardScalerを利用すれば簡単である。

　分析データはsklearnにあるboston近郊の住宅価格に関するデータを使用した。説明変数Xは13因子ある。まず、主成分分析の結果を下図に示す。縦軸は全データの分散に対する寄与率で、横軸は左から順に、第一主成分、第二主成分、…の主成分を表し、PCはPrincipal Componentの略である。データをみると、第一主成分PC1で47％を占める。桃色の折れ線は累積寄与率であり80％を超えるのはPC5である。この段階で次元削減は厳しいとは思われる。PC2～3で80％越えになるような取り扱いデータの場合に効果的なのではと思う。

そして、下図に横軸にPC1、縦軸にPC2をとり、全データをプロットしたグラフを示す。赤矢印は各説明変数XのPC1とPC2に対する方向性を示し、その長さは影響度の高さを示す。グラフの見方だが、横軸PC1の意味は全データ中の最大の分散方向を表すため、PC1軸方向にベクトルが長いものが影響度が高いととれる。ここで、注意点は全データの生点が密集している方向に指してる説明変数Xほど、精度が高いととらえる。ベクトルが点の少ない方向を示している説明変数Xの場合は、データ点数が少ないことに依る飛び値の影響を受けている可能性が考えられる。下図から読み取れることは、左上を向いているBとRM、右下を向いてるLSTAT、右上を向いてるNOXは説明変数として少なくとも必要だろうと考える。更に、BとRMはベクトル方向がほぼ同じなため、両者には依存性があるのではないか等と考えたりする。

上記のような考え方から、説明変数13因子中から8因子を除いた5因子の行列散布図を下図に示す。解釈が妥当かは絞り込んだ散布図でも確認することは必要に思う。

そして、説明変数Xを13因子から5因子に絞った条件で重回帰分析した結果を下図に示す。青点は生データであり、赤点は得られた回帰式に生データXを代入したものである。

数値でみると、決定係数は下図左の74％から下図右の67％へと低下した。上図右上の箇所が青点と赤点の重なりが少なく、回帰式で再現性が低い様子がうかがえる。

　

一方、13因子で重回帰分析したグラフ結果を下図に示す。図中の右上箇所の青点と赤点の重なりは先のグラフと比較して多い。

▼本プログラム

データの標準化、主成分分析、重回帰分析のいずれもコードは2～3行だが、視覚化するためのグラフ化の方が凡そ5～10行でやや長い。

# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import seaborn as sns

from sklearn import datasets
from sklearn import linear_model
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

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
    plt.figure(figsize=(10,10))
    sns.heatmap(x_sc.corr(), annot=True)
    plt.show()
    return x_sc

def my_pca(X_sc, Y): # 主成分分析
    ################################# Explained Variance Rate
    pca = PCA(n_components=13)
    pca.fit(X_sc)
    ### plot
    x = ['PC%02s' %i for i in range(1, len(pca.explained_variance_ratio_)+1)]
    y = pca.explained_variance_ratio_
    cum_y = np.cumsum(y)
    plt.figure(figsize=(9,4))
    plt.bar(x, y, align="center", color="green")
    plt.plot(x, cum_y, color="magenta", marker="o")
    for i, j in zip(x, y):
        plt.text(i, j, '%.2f' % j, ha='center', va='bottom', fontsize=14)
    plt.ylim([0,1])
    plt.ylabel('Explained Variance Rate', fontsize = 14)
    plt.tick_params(labelsize = 14)
    plt.tight_layout()
    plt.grid()
    plt.show()
    ################################# PC1 vs PC2 (PC:Principal Component)
    x_pc = pca.fit_transform(X_sc)
    ### plot
    x_data = x_pc[:,0]
    y_data = x_pc[:,1]
    pc0 = pca.components_[0]
    pc1 = pca.components_[1]
    plt.figure(figsize=(7,7))
    plt.scatter(x_data, y_data, c=Y/len(set(Y)), marker=".")
    for i in range(pc0.shape[0]):
        plt.arrow(0, 0, pc0[i]*8, pc1[i]*8, color='r', width=0.0005, head_width=0.1, alpha=0.9)
        plt.text(pc0[i]*9, pc1[i]*9, X_sc.columns.values[i], color='r')
    plt.xlim([-6,6])
    plt.ylim([-6,6])
    plt.xlabel('Principal Component 1', fontsize = 14)
    plt.ylabel('Principal Component 2', fontsize = 14)
    plt.tick_params(labelsize = 14)
    plt.grid()
    c = patches.Circle(xy=(0, 0), radius=2, ec='r', fill=False)
    plt.axes().add_patch(c)
    plt.show()

def my_regression(X, Y): # 線形重回帰分析
    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('raw X-data')
    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='Multiple regression analysis')
    plt.legend(loc='lower right', fontsize=12)
    plt.show()

if __name__ == '__main__':
    # データセットを読み込みXとYで返す
    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)
    
    # 主成分分析
    my_pca(X_sc, Y)
    
    # 説明変数Xを減らす
    # 　主成分分析から読み取って独自判断で説明変数Xを削除
    drop_col_1 = ['ZN', 'PTRATIO', 'CRIM', 'RAD', 'INDUS', 'AGE', 'CHAS', 'B']
    X1 = X.drop(drop_col_1, axis=1)        
    # 　全データ用いての回帰結果の係数から感度の低い説明変数Xを削除
    drop_col_2 = ['LSTAT', 'CRIM', 'TAX', 'AGE', 'B', 'INDUS', 'ZN', 'RAD']
    X2 = X.drop(drop_col_2, axis=1)
    ### Matrix scatter plot
    df_X1Y = pd.concat([X1,Y], axis=1)
    sns.pairplot(df_X1Y)
    plt.show()
    df_X2Y = pd.concat([X2,Y], axis=1)
    sns.pairplot(df_X2Y)
    plt.show()

    # 線形重回帰分析
    my_regression(X, Y)  # 1.全データを使用した場合の回帰
    my_regression(X1, Y) # 2.主成分分析から読み取って独自判断で説明変数Xを減らしての回帰
    my_regression(X2, Y) # 3.1の回帰結果の係数から感度の低い説明変数Xを減らしての回帰

　主成分分析が効果を発揮するのは、扱う多変量データの説明変数X間で分散の差が大きく、且つそれらが出来るだけ連続的に多数のデータが分布している場合であると思う。データ数が少なく飛び値的な理由で分散が大きいデータが含まれている場合には判断を見誤る可能性がある。これを防ぐためには散布図も併用して確認する必要がある。例えばこのbostonデータセットにはダミー変数が設定されている。CHASの項は川に近い場合は1、遠い場合は0であり不連続データであるため、飛び値同様の症状が起こりうる。

 結局は扱うデータに依るため、実施してみないとわからない。善し悪しの判断は、分類や回帰をしてからの使用目的に依ると考える。

参考リンク

sklearn.decomposition.PCA — scikit-learn 0.19.1 documentation

sklearn.preprocessing.StandardScaler — scikit-learn 0.19.1 documentation

patches — Matplotlib 2.2.2 documentation

hk29.hatenablog.jp

hk29.hatenablog.jp

以上

<広告>

リンク