本家ブログ（実践ケモインフォマティクス）もよろしくお願いします。

主成分分析について

概要

教師なし手法の一つであり、PCA（Principal Component Analysis）と呼ばれる。
多次元のデータを低次元化する手法であり、データの可視化によく用いられる。（低次元化し2次元にマッピングする事でデータの可視化が可能）
データセットの情報量損失を最小限にしつつ低次元化することが可能。
データのノイズを寄与率の小さい主成分として除くことが可能。
PCA後の変数（主成分）は互いに相関0なので、これを用いて回帰分析すると回帰係数の値が安定になる。（この手法は主成分回帰（PCR）と呼ばれる）

方法（ざっくり）

データのばらつきが最も大きくなる方向に新しく軸を設定する。
1で設定した軸と直行し、かつデータのばらつきが最も大きくなる方向に2軸目を設定する。
1,2を繰り返す。
1,2の過程では、それぞれの主成分がデータ全体をどの程度表せているかの指標（寄与率）が得られる。

主成分分析で得られる指標

固有値：各主成分に対応した固有値であり、主成分の分散に対応する。
寄与率：ある主成分の固有値（分散）が表す情報が、データの全ての情報（分散）の何割を表しているかの指標。（言い換えると、ある主成分がデータの何割を表せられているかの指標のこと）
累積寄与率：各主成分の寄与率を大きい順に足したもの。そこまでの主成分でデータの何割を表現できているかを表す。（累積寄与率が70〜80％がカバーできていればOKと判断する場合が多い）

注意点

主成分分析の前に変数を標準化する必要がある。

コード

データの準備

水溶解度のデータセットを利用する（金子研究HP から取得可能）

import os
import pandas as pd
import numpy as np
import math
import matplotlib.figure as figure
import matplotlib.pyplot as plt
from rdkit import Chem
from rdkit.ML.Descriptors import MoleculeDescriptors

# データの読み込み
suppl = Chem.SDMolSupplier('logSdataset1290_2d.sdf')
mols = [mol for mol in suppl if mol is not None]
# RDkit記述子の計算

descriptor_names = [descriptor_name[0] for descriptor_name in Chem.Descriptors._descList]
descriptor_calculation = MoleculeDescriptors.MolecularDescriptorCalculator(descriptor_names)
RDkit = [descriptor_calculation.CalcDescriptors(mol) for mol in mols]
# データフレーム作成

df = pd.DataFrame(RDkit, columns = descriptor_names)
df['logS'] =  [mol.GetProp('logS') if 'logS' in mol.GetPropNames() else 'None' for mol in mols]
df = df.astype(float)

# サンプル数を少なくする。
df.sort_values('logS',ascending=False,inplace=True)
df = df.iloc[::10,:]

# 相関係数が高い10変数を選抜する。
df_corr = df.corr()
df_corr.sort_values('logS',ascending=False,inplace=True)
selected_features = df_corr.index[0:11]
print(selected_features)
df = df.loc[:,selected_features]

df_X = df.drop('logS',axis=1)
df_y = df['logS']

主成分分析の実施とデータの可視化

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA #主成分分析器
import matplotlib.ticker as ticker

# 標準化
scaler_X = StandardScaler()
X_scaled = scaler_X.fit_transform(df_X)

# 主成分分析の実行
pca = PCA()
feature = pca.fit_transform(X_scaled)
df_pca = pd.DataFrame(feature, index=df_X.index,
columns=["PC{}".format(x+1) for x in range(feature.shape[1])])

# データを主成分空間に写像
plt.figure(figsize=(5,4))
plt.scatter(df_pca.iloc[:, 0],df_pca.iloc[:, 1], alpha=0.8,
c=df_y, vmin=df_y.min(),vmax=df_y.max())
plt.colorbar()

# 図にテキストを表示する場合
# for x, y, name in zip(df_pca.iloc[:, 0], df_pca.iloc[:, 1], df_pca.index.tolist()):
#     plt.text(x, y, name)
plt.grid()
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.show()

# 累積寄与率を図示する
plt.figure(figsize=(4,4))
plt.gca().get_xaxis().set_major_locator(ticker.MaxNLocator(integer=True))
plt.plot([0] + list(np.cumsum(pca.explained_variance_ratio_)), "-o")
plt.xlabel("Number of principal components")
plt.ylabel("Cumulative contribution rate")
plt.grid()
plt.show()

第1,第2主成分でのマッピング結果

累積寄与率の表示

詳細な計算式など

計算式（ざっくり）

以下数式は金子研究室HP記載内容を参考にしています。より詳細な計算は金子研究室HPやqiitaの記事についてがおすすめです。
また、数式に間違いがあればコメントくだされば助かります。

主成分の計算

以下のようなデータを使った、2変数のPCAを考える.

	変数A	変数B
サンプル1	$x_A^{(1)}$	$x_B^{(1)}$
サンプル2	$x_A^{(2)}$	$x_B^{(2)}$
$\vdots$	$\vdots$	$\vdots$
サンプルn	$x_A^{(n)}$	$x_B^{(n)}$

ここで、

$$ \boldsymbol{X_A} = \begin{bmatrix} x_A^{(1)} \\ x_A^{(2)} \\ \vdots \\ x_A^{(n)} \end{bmatrix} \qquad \boldsymbol{X_B} = \begin{bmatrix} x_B^{(1)} \\ x_B^{(2)} \\ \vdots \\ x_B^{(n)} \end{bmatrix} $$

と表記する。
第 $i$ 主成分を $t_i$、第 $i$ 主成分に対応する$j$番目の変数の重み（ローディング）を$p_i^{(j)}$とすると、

$$ t_1 = \boldsymbol{X}_Ap_1^{(1)} + \boldsymbol{X}_Bp_1^{(2)} \\ t_2 = \boldsymbol{X}_Ap_2^{(1)} + \boldsymbol{X}_Bp_2^{(2)}$$

これを行列で表すと、

$$ \begin{bmatrix} t_1^{(1)} & t_2^{(1)}\\ t_1^{(2)} & t_2^{(2)}\\ \vdots & \vdots\\ t_1^{(n)} & t_2^{(n)} \end{bmatrix}= \begin{bmatrix} x_A^{(1)} & x_B^{(1)} \\ x_A^{(2)} & x_B^{(2)} \\ \vdots & \vdots \\ x_A^{(n)} & x_B^{(n)} \end{bmatrix} \begin{bmatrix} p_1^{(1)} & p_2^{(1)} \\ p_1^{(2)} & p_2^{(2)} \\ \end{bmatrix}$$

$x_i^{(k)}$ ：サンプル k における、 i 番目の変数の値
$t_i^{(k)}$ ：サンプル k における、第 i 主成分の値
$p_i^{(j)}$ ：第 i 主成分に対応する、 j 番目の変数の重み

主成分分析では、データセットのばらつき（分散）が最大の方向を主成分軸とする。
ここで、元のデータの平均が0であれば、軸をとりなおした後のデータの平均も0となることを利用すると、
軸変更後のデータの分散が最大となるためには、主成分の値の二乗和が最大になれば良い。
（第1主成分について考えると）つまり、

$$ S = \sum_{i=1}^{n}(t_1^{(i)})^2 $$

ここで、ローディング（重み）によって、主成分が変わらないように、以下の制約を設ける。

$$ (p_1^{(1)})^2 + (p_1^{(2)})^2 = 1　$$

上記の制約と、Lagrangeの未定乗数法を利用して $S$ が最大となる条件を求めると、

$$ (\boldsymbol{X}^T\boldsymbol{X} - \lambda\boldsymbol{E})\boldsymbol{p_1} = 0　$$

$$ \boldsymbol{X} = \begin{bmatrix} x_A^{(1)} \\ x_A^{(2)} \\ ⋮\\ x_A^{(n)} \end{bmatrix} , \quad \boldsymbol{E} = \begin{bmatrix} 1&0\\ 0&1\\ \end{bmatrix} , \quad \boldsymbol{p_1} = \begin{bmatrix} p_1^{(1)}\\ p_1^{(2)}\\ \end{bmatrix} \qquad　$$

となる。
この式が $p_1^{(1)}=p_1^{(2)}=0$ 以外の解を持つためには、$(\boldsymbol{X}^T\boldsymbol{X} - \lambda\boldsymbol{E})$　の行列式が０である必要がある。

つまりこの問題は、$\lambda$ を固有値、$\boldsymbol{p_1}$　及び、$\boldsymbol{p_2}$ を固有ベクトルとする固有値問題と考えることができる。
これによって $\boldsymbol{p_1}$　、$\boldsymbol{p_2}$ を求め、対応する主成分を計算する。

(累積)寄与率の計算

第 $i$ 主成分に対応する固有値 $\lambda_i$ は、その主成分の二乗和（=分散）に等しいので、固有値 $\lambda_i$ を第 $i$ 主成分の持つ情報量と仮定すると、寄与率 $c_i$ は以下の式で求められる。

$$ c_i = \frac{\lambda_i}{\sum_{j=1}^m\lambda_j}　$$

（$m$ ：全ての主成分の数）

また、第 $i$ 主成分までの寄与率の和は、第 $i$ 主成分までの累積寄与率 $C_i$ と呼ばれる。

$$ C_i = \frac{\sum_{j=1}^i\lambda_j}{\sum_{j=1}^m\lambda_j}　$$

計算コード

2変数に絞って、sklearnを使わずに主成分分析してみる。

データの準備

# データの抜き出しと標準化
df_X2 = df_X.iloc[:,:2]
scaler_X = StandardScaler()
X_scaled2 = scaler_X.fit_transform(df_X2)

# 変換前のデータを可視化する。
plt.figure(figsize=(5,5))
plt.scatter(X_scaled2[:,0],X_scaled2[:,1], alpha=0.8)
plt.xlabel("X1")
plt.ylabel("X2")
plt.grid()
plt.show()

主成分分析前のデータはこんな感じ

主成分分析実行

$\boldsymbol{X}^T\boldsymbol{X}$ の固有ベクトル(Eigenvector）と固有値(eigenvalue)を求める。

import numpy as np
import numpy.linalg as LA

# 固有値、固有ベクトルの計算
a= np.dot(X_scaled2.T, X_scaled2)
eigenvalue, eigenvector = LA.eig(a)
print('固有値:\n',eigenvalue)
print('----')
print('固有ベクトル:\n',eigenvector)
"""output
固有値:
 [246.22143135  11.77856865]
----
固有ベクトル:
 [[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
"""

# 主成分の計算
pca = np.dot(X_scaled2, eigenvector)

# 寄与率の計算
eigenvalue_sum = eigenvalue.sum()
eigenvalue1 = eigenvalue[0]/eigenvalue_sum
eigenvalue2 = eigenvalue[1]/eigenvalue_sum
print('----')
print('寄与率(第1主成分):',eigenvalue1)
print('寄与率(第2主成分):',eigenvalue2)
"""output
寄与率(第1主成分): 0.9543466331268278
寄与率(第2主成分): 0.04565336687317221
"""

# 変換のデータ
plt.figure(figsize=(5,5))
plt.scatter(pca[:,0],pca[:,1], alpha=0.8)
plt.xlabel("pc1")
plt.ylabel("pc2")
plt.grid()
plt.savefig('after_pca.png')
plt.show()

主成分分析後のマッピング

参考