はじめに
scikit-learn（sklearn）での主成分分析（PCA）の実装について解説していきます。
- Pythonで主成分分析を実行したい方
- sklearnの主成分分析で何をしているのか理解したい方
- 主成分分析の基本中の基本（.fitや.transform）プラスアルファを学びたい方
の参考になれば嬉しいです。この投稿ではPythonでの実装に焦点を当てますので、理論的な内容を学びたいという方は是非こちらの投稿を参考にして下さい。
Pythonでの実装
前半はデータの前処理について説明し、後半から主成分分析について解説していきます。「とにかく主成分分析自体の実装方法が知りたい！」という方は主成分分析の実行までスキップして下さい。
環境
- conda 4.8.3
- python 3.8.3
- pandas 1.0.5
- matplotlib 3.2.2
- numpy 1.18.5
- scikit-learn 0.23.1
データの前処理
主成分分析ができるように説明変数を定量データに絞り、不要な特徴量を削除していきます。データはこちらに保存しているcompound_2.csvを使用します。コードは「化学のための Pythonによるデータ解析・機械学習入門」を参考に書かせていただいている部分が多いです。
データの読み込み、定性データの削除
データを読み込みます。compound_2.csvのファイルは実行中のディレクトリに保存する必要があります。
import pandas as pd
dataset = pd.read_csv('compound_2.csv', index_col=0)
dataset.head()
- index_colでindexとして使いたい列を指定します
以降、**band_gapを目的変数として扱っていきます。**読み込んだデータセットを編集していきます。band_gap＞0のデータに限定、インデックス名の変更、定性データの削除を行います。
# band_gapが正のデータに限定する
dataset = dataset[dataset["band_gap"] > 0]
# 化合物名（pretty_formula列）をindexに採用する（listとして引き渡さないと行として認識されるので注意）
dataset.index = list(dataset.loc[:,"pretty_formula"])
# データを定量データのみにする
dataset = dataset[["band_gap", "density", "e_above_hull", "energy", "energy_per_atom", "formation_energy_per_atom", "nelements", \
"nsites", "total_magnetization", "volume"]]
dataset.head()
だいぶスッキリしてきましたね。ここで、現状のデータの大きさを確認しておきます。
# データセットの大きさの確認
dataset.shape
(229, 10)
特徴量の確認（相関行列、散布図）
相関行列や散布図を用いて特徴量の分布などを確認していきます。
# ライブラリのインポート
import matplotlib.pyplot as plt
import seaborn as sns
# 相関行列
correlation_coefficients = dataset.corr() # 相関行列の計算
# 相関行列のヒートマップ (相関係数の値あり)
plt.rcParams['font.size'] = 12
plt.figure(figsize=(12, 8)) # この段階で画像のサイズを指定する
sns.heatmap(correlation_coefficients, vmax=1, vmin=-1, cmap='seismic', square=True, annot=True, xticklabels=1, yticklabels=1)
plt.xlim([0, correlation_coefficients.shape[0]])
plt.show()
相関係数のヒートマップより、volumeとnsitesは強い正の相関がありenergyとnsites,volumeは強い負の相関があることが読み取れます。nsitesの物理的な解釈ですが、こちらを参考にするとlattice position(結晶格子点)の数ではないかと考えられます。説明変数の間に強い相関がある際は多重共線性によりオーバーフィットする可能性が高くなるので、volumeと強い正の相関があるnsitesを変数として削ることにします。
続いて、特徴量の散布図を描いていきます。
plt.rcParams['font.size'] = 8 # 横軸や縦軸の名前の文字などのフォントのサイズ
pd.plotting.scatter_matrix(dataset, c='blue', figsize=(10, 10)) # ここで画像サイズを指定
plt.show()
「volumeとnsitesは強い正の相関がある」ことは述べましたが、散布図（一番下の列の右から3番目）を見てもこの傾向が視覚的に読み取れます。nelementsは一定の値しかとらない（分散が0である）ので変数として削ることにします。こちらにも書かれていますが、分散が0の変数（一定値しかとらない変数）は特徴量として意味をなさないので削除します。
# nelementsとnsitesを説明変数から削る
dataset = dataset.drop(['nelements', 'nsites'], axis=1) # axis=1で横方向に見ることを指定して特定列を削除
dataset.head()
このように機械学習を適用する際は特徴量をなるべく減らすということがモデルの精度向上や計算量削減の観点から定石となっております。今回は以下の2つを実行して特徴量を減らしました。
- 2変数間に強い相関がある際はどちらか片方に絞る
- 分散が0（一定値しかとらない）の変数は削除する
ここまで加工したデータをcsvファイルに書き出して保存しておきます。
dataset.to_csv("compound_2_edit.csv", encoding="shift_jis")
主成分分析の実行
データの読み込み
上記で加工したデータを読み込み、データの大きさを確認します。なお、加工済みのデータはこちらにcompound_2_edit.csvとして保存しています。
dataset = pd.read_csv('compound_2_edit.csv', index_col=0)
dataset.shape
(229, 8)
主成分分析を実行する前にデータを標準化します。今回のデータでは1列目のband_gapを目的変数とし、2列目のdensity以降を説明変数として扱っていきます。
x = dataset.iloc[:, 1:] # 最初の列のbandgapを目的変数とする
autoscaled_x = (x - x.mean()) / x.std() # オートスケーリング
autoscaled_x
モデルの宣言
主成分分析を実行していきます。公式ページも分かりやすく、参考になります。
from sklearn.decomposition import PCA
pca = PCA() # PCA を行ったり PCA の結果を格納したりするための変数を、pca として宣言
pca.fit(autoscaled_x) # PCA を実行
PCA()
PCA()のパラメータとして一般的なのは"n_components"であり、主成分数を定義します。何も指定しない際は全ての成分数が保持されます。（つまり、今回であれば主成分数は7として扱われます）
ローディング（各主成分の構成方法）
**ローディングとは各主成分がどのように構成されるかを表す値です。**以下のコードでローディングを確認します。
loadings = pd.DataFrame(pca.components_.T, index=x.columns)
loadings
第一主成分は(-0.24)×(density) + (0.06)×(e_above_hull) + … + (0.68)×(volume) で構成されていることが読み取れます。今回はn_componentsを指定していないので、もともとの特徴量の数と同じ7個の主成分が求められます。
スコア（各主成分に対するサンプルの値）
**スコアとは各サンプルが各主成分軸上のどの座標に位置するかを表す値です。**以下のコードでスコアを確認します。
score = pd.DataFrame(pca.transform(autoscaled_x), index=x.index)
score
LiHは第一主成分軸上では-0.81、第二主成分軸上では-2.55に位置すると読み取れます。
ここでデータフレームの最初の4つのデータ（LiH, BeH2, B9H11, B4H5）が第一主成分vs第二主成分平面においてどこに位置するかを可視化してみます。
num = 4 # 可視化するデータ数を指定
plt.scatter(score.iloc[:num,0], score.iloc[:num,1])
plt.rcParams["font.size"] = 11
# プロットしたデータにサンプル名をラベリング
for i in range(num):
plt.text(score.iloc[i,0], score.iloc[i,1], score.index[i], horizontalalignment="center", verticalalignment="bottom")
plt.xlim(-1, 2.5)
plt.ylim(-3, -1)
plt.yticks(np.arange(-3, -0.5, 0.5))
plt.xlabel("t1")
plt.ylabel("t2")
plt.grid()
plt.show()
横軸は第一主成分を表し、縦軸は第二主成分を表します。確かにLiHは第一主成分軸上では-0.81、第二主成分軸上では-2.55に位置していることが分かります。
寄与率（表現できる情報量の割合）
**寄与率とは表現できる情報量の割合を表す値です。**各成分の寄与率と累積寄与率を求めていきます。
# 寄与率を算出
contribution_ratios = pd.DataFrame(pca.explained_variance_ratio_)
contribution_ratios
第一主成分だけでもともとの特徴量全体で表していた情報（7つの説明変数を使って表現していた情報）の29%を表現できていることが分かります。
# 累積寄与率を算出（.cusum()で累積和を計算 .sum()では総和しか得られない）
cumulative_contribution_ratios = contribution_ratios.cumsum()
cumulative_contribution_ratios
第一から第四主成分までの4つの主成分を用いることでもともとの情報量の85%が表現できます。また、もとの特徴量と同じ数の主成分を使えば情報量の損失はないので主成分数7で累積寄与率は1となります。主成分数に対する寄与率を以下のコードで可視化します。
cont_cumcont_ratios = pd.concat([contribution_ratios, cumulative_contribution_ratios], axis=1).T
cont_cumcont_ratios.index = ['contribution_ratio', 'cumulative_contribution_ratio'] # 行の名前を変更
# 寄与率を棒グラフで、累積寄与率を線で入れたプロット図を重ねて描画
x_axis = range(1, contribution_ratios.shape[0] + 1) # 1 から成分数までの整数が x 軸の値
plt.rcParams['font.size'] = 18
plt.bar(x_axis, contribution_ratios.iloc[:, 0], align='center') # 寄与率の棒グラフ
plt.plot(x_axis, cumulative_contribution_ratios.iloc[:, 0], 'r.-') # 累積寄与率の線を入れたプロット図
plt.xlabel('Number of principal components') # 横軸の名前
plt.ylabel('Contribution ratio(blue),\nCumulative contribution ratio(red)') # 縦軸の名前。\n で改行しています
plt.show()
データの可視化
主成分分析の結果を可視化します。第一主成分vs第二主成分平面に各サンプルをプロットします。また、各プロットはband_gap（目的変数）の値で色付けしておきます。
# 第 1 主成分と第 2 主成分の散布図 (band_gap の値でサンプルに色付け)
plt.scatter(score.iloc[:, 0], score.iloc[:, 1], c=dataset.iloc[:, 0], cmap=plt.get_cmap('jet'))
clb = plt.colorbar()
clb.set_label('band_gap', labelpad=-20, y=1.1, rotation=0)
plt.xlabel('t1')
plt.ylabel('t2')
plt.show()
この散布図より、第一主成分の値が大きいとband_gapが小さくなる傾向が読み取れます。また、プロットが黄色や赤色のband_gapの大きい材料は第一主成分が0以下、第二主成分が2以下の領域に多いことも読み取れます。なお、カラーバーの表示方法はこちらを参考にしました。
このように主成分分析によって次元削減やデータの可視化ができます。（実際にはこの後の仕事として、ローディングの値から目的変数に影響を与える因子を抽出したり、主成分分析で生成した主成分を使って機械学習モデルの作成を行ったりすると考えられます）
まとめ
scikit-learn（sklearn）での主成分分析（PCA）の実装について解説致しました。（少し饒舌になり過ぎたでしょうか…）