主成分分析の概要とpython実装

[本家ブログ（実践ケモインフォマティクス）](https://datadriven-rnd.com/2021-02-08-222511/) もよろしくお願いします。

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


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

![Fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/348376/1566613e-4004-db8c-1716-737274bb7368.png)

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


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

## コード
### データの準備

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

```python
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']
```
### 主成分分析の実施とデータの可視化

```python
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主成分でのマッピング結果**
![PCA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/348376/f1bb1bd3-d8c1-d8a1-b1d4-15f34e486d13.png)

**累積寄与率の表示**
![Cumulative_contribution_rate2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/348376/118c0ff8-b71a-a8e0-27e8-44600b9f2ab4.png)

## 詳細な計算式など
### 計算式（ざっくり）

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

### 主成分の計算

以下のようなデータを使った、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)}$ |

ここで、

<p>
$$ \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} $$
</p>

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

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

これを行列で表すと、

<p>
$$ \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}$$
</p>

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


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

<p>$$ S = \sum_{i=1}^{n}(t_1^{(i)})^2 $$</p>

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

<p>$$ (p_1^{(1)})^2 + (p_1^{(2)})^2 = 1　$$</p>

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

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

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

となる。
この式が $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$ は以下の式で求められる。

<p>$$ c_i = \frac{\lambda_i}{\sum_{j=1}^m\lambda_j}　$$</p>
（$m$ ：全ての主成分の数）

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

<p>$$ C_i = \frac{\sum_{j=1}^i\lambda_j}{\sum_{j=1}^m\lambda_j}　$$</p>

## 計算コード
2変数に絞って、sklearnを使わずに主成分分析してみる。

#### データの準備
```python
#データの抜き出しと標準化
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()
```

主成分分析前のデータはこんな感じ
![before_pca.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/348376/2333f5e6-85d3-4dfd-abf9-9ed782cdbb97.png)

#### 主成分分析実行
$\boldsymbol{X}^T\boldsymbol{X}$ の固有ベクトル(Eigenvector）と固有値(eigenvalue)を求める。

```python
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()
```

主成分分析後のマッピング
![after_pca.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/348376/0ec22a0c-4ba3-fe6d-7d42-0e2ec056f5b0.png)


## 参考
- [主成分分析(Principal Component Analysis, PCA)～データセットの見える化・可視化といったらまずはこれ！](https://datachemeng.com/principalcomponentanalysis/)
- [主成分分析(Principal Component Analysis, PCA)の前に変数の標準化(オートスケーリング)をしたほうがよいのか？](https://datachemeng.com/autoscaling_before_pca/)
- [主成分分析](https://www.macromill.com/service/data_analysis/principal-component-analysis.html)
- [主成分分析のまとめ](https://qiita.com/amateur2020/items/59594fa66c4a46e61906)