化合物の溶解度予測（予測の信頼性も評価する）

# はじめに
[化合物の溶解度予測（基礎編）](https://qiita.com/ottantachinque/items/f172bdf771ac87f1a962)をベースにして、今回は**予測の信頼性（=予測値の分散=予測値がどれだけ不安定か）**も合わせて評価する。


# 目的
分子の水溶解度（LogS）の予測および、予測値の分散を求める。

#諸条件

##データ

100分子の水溶解度（LogS）データ（sdfファイル形式）
（[金子研究室HPのデータセット](https://datachemeng.com/pythonassignment/)から取得した1290分子のデータの中からランダムにピックアップ）

##変数作成
- RDkit 記述子を使用（[参考HP](https://future-chem.com/rdkit-descriptors/#RDKit)）

##前処理
- 以下の変数を削除
  - 分散が0の変数
  - 互いに相関係数が高い（0.95以上）変数ペアの内、片方の変数  
    （金子研モジュール使用）
  - 50%以上の分子(50分子）で同じ値をとる変数  
  (テストデータの割合を50%に設定しているため） 
- 変数は全て標準化（平均：0, 分散：1）

##解析手法

- 予測器：ランダムフォレスト
- ハイパーパラメータ最適化：なし（省略のため）

###予測値の平均・分散の計算方法
以下の手順(1,2)を200回繰り返し、各分子の予測値の平均値、標準偏差を求めた。
（散布図では ±2σ の範囲をエラ〜バーで表示した）
 1. 訓練データとテストデータをランダムに分割（テストデータ比率50%）
 2. 訓練データで予測器を学習させ、テストデータでの予測値を求める。

#前準備

- データの取得（[金子研究室HPのデータセット](https://datachemeng.com/pythonassignment/)を利用）
- モジュールのインストール

```
conda install -c conda-forge rdkit
pip install dcekit
``` 

#コード

###モジュールのインポート

```python
import os
import pandas as pd
import numpy as np
import math

from scipy.spatial.distance import cdist

import matplotlib.figure as figure
import matplotlib.pyplot as plt

from rdkit import Chem
from rdkit.ML.Descriptors import MoleculeDescriptors

from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor as rf
from sklearn.model_selection import train_test_split,GridSearchCV
from sklearn.metrics import r2_score, mean_absolute_error,mean_squared_error

from dcekit.variable_selection import search_high_rate_of_same_values, search_highly_correlated_variables
```

###データの読み込みと変数作成


```python
#データの読み込み
suppl = Chem.SDMolSupplier('logSdataset1290_2d.sdf')
mols = [mol for mol in suppl if mol is not None]

#SDFファイル中のプロパティの取得
property_list = list(mols[0].GetPropNames())
print(property_list)
# =>['CAS_Number', 'logS']

#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)
#溶解度(logS)の値を取得
df['logS'] =  [mol.GetProp('logS') if 'logS' in mol.GetPropNames() else 'None' for mol in mols]
```

###変数のスクリーニング

```python
df = df.astype(float)

#ランダムに100サンプルを抽出
df = df.sample(n=100, random_state=0)

#欠損のある変数を削除
df.dropna(axis=1,inplace=True)

#分散が０の変数削除
del_num1 = np.where(df.var() == 0)
df = df.drop(df.columns[del_num1], axis=1)

#変数選択（互いに相関関係にある変数の内一方を削除）   
threshold_of_r = 0.95 #変数選択するときの相関係数の絶対値の閾値
corr_var = search_highly_correlated_variables(df, threshold_of_r)
df.drop(df.columns[corr_var], axis=1, inplace=True)

#同じ値を多くもつ変数の削除
inner_fold_number = 2 #CVでの分割数（予定）
rate_of_same_value = []
threshold_of_rate_of_same_value = (inner_fold_number-1)/inner_fold_number
for col in df.columns:
    same_value_number = df[col].value_counts()
    rate_of_same_value.append(float(same_value_number[same_value_number.index[0]] / df.shape[0]))
del_var_num = np.where(np.array(rate_of_same_value) >= threshold_of_rate_of_same_value)
df.drop(df.columns[del_var_num], axis=1, inplace=True)

print(df.shape)
# => (100, 53)
```

###学習と予測
訓練・テストデータの分割方法を変えて、学習と予測を繰り返す

```python
df_X = df.drop(['logS'],axis=1)
df_y = df.loc[:,['logS']]

#予測値を保存するDataFrame
df_predict = pd.DataFrame(index=df.index)

for i in range(200):

　　#訓練、テストデータの作成
    df_X_train, df_X_test, df_y_train, df_y_test = train_test_split(df_X, df_y, test_size=0.5, random_state=i)
    
    X_train = df_X_train.values
    y_train = df_y_train.values
    X_test = df_X_test.values
    y_test = df_y_test.values
    
　　#標準化
    scaler_X = StandardScaler()
    scaler_y = StandardScaler()

    scaled_X_train = scaler_X.fit_transform(X_train)
    scaled_y_train = scaler_y.fit_transform(y_train)
    scaled_X_test = scaler_X.transform(X_test)
    scaled_y_test = scaler_y.transform(y_test)

　　#学習と予測
    model = rf()
    model.fit(scaled_X_train, scaled_y_train)

    scaled_pred_y_test = np.ndarray.flatten(model.predict(scaled_X_test))
    pred_y_test = scaler_y.inverse_transform(scaled_pred_y_test)
    
    df_predict.loc[df_X_test.index, 'pred_y_{}'.format(i)] = pred_y_test
```

###予測結果の可視化

```python
#予測値の平均と分散を計算
pred_y_test = df_predict.mean(axis=1).values
pred_y_test_std = df_predict.std(axis=1).values * 2 
y_test = df_y.values

#精度評価指標を計算
print('r2_test:{}'.format(round(r2_score(y_test, pred_y_test),3)))
print('MAE_test:{}'.format(round(mean_absolute_error(y_test, pred_y_test),3)))
print('RMSE_test:{}'.format(round(np.sqrt(mean_squared_error(y_test, pred_y_test)),3)))
"""
(output)
r2_test:0.819
MAE_test:0.667
RMSE_test:0.847
"""

#結果の可視化（yyplot作成）
plt.figure(figsize=(6,6))
plt.errorbar(y_test, pred_y_test, yerr=pred_y_test_std, capsize=5, alpha=0.5,
             fmt='o', 
             markersize=10, 
             ecolor='black', 
            )
y_max_ = max(y_test.max(), pred_y_test.max())
y_min_ = min(y_test.min(), pred_y_test.min())
y_max = y_max_ + (y_max_ - y_min_) * 0.1
y_min = y_min_ - (y_max_ - y_min_) * 0.1

plt.plot([y_min , y_max],[y_min, y_max], 'k-')

plt.ylim(y_min, y_max)
plt.xlim(y_min, y_max)
plt.xlabel('logS(observed)',fontsize=20)
plt.ylabel('logS(predicted)',fontsize=20)
plt.title('yyplot(logS)',fontsize=20)

plt.show()
```
![yyplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/348376/8edca0c3-8edd-1adf-5558-8fb4e2adf12a.png)



##サンプル同士の距離と予測値のばらつきの関係
[k-最近傍法](https://datachemeng.com/knn/)を使って、各サンプルと他の全サンプルの距離の計算を計算して、予測値のばらつきとの関係を調べた。


```python
#距離計算用の変数セット（標準化したもの）を用意
scaler_X_all = StandardScaler()
df_X_scaled = pd.DataFrame(scaler_X_all.fit_transform(df_X), index=df_X.index, columns=df_X.columns)

#各サンプル同士の距離行列作成し、全てのサンプルとの距離の平均を計算する
distance_between_autoscaled_df = cdist(df_X_scaled, df_X_scaled)
knn_distance = np.mean(distance_between_autoscaled_df, axis=1)

#平均距離の情報と予測値の標準偏差をpandas DataFrame に格納
df_distance = pd.DataFrame({'distance':knn_distance, 'std':df_predict.std(axis=1).values}, index=df.index)
df_distance['base'] = 0

#以下は結果の可視化
plt.figure(figsize=(6,6))
plt.errorbar(df_distance['distance'],df_distance['base'], yerr=df_distance['std'], capsize=5, alpha=0.5,
             fmt='o', 
             markersize=10, 
             ecolor='black', 
            )

y_max_ = df_distance['std'].max()
y_min_ =  - df_distance['std'].max() 
X_max_ = df_distance['distance'].max()
X_min_ = df_distance['distance'].min()

y_max = y_max_ + (y_max_ - y_min_) * 0.1
y_min = y_min_ - (y_max_ - y_min_) * 0.1
X_max = X_max_ + (X_max_ - X_min_) * 0.1
X_min = X_min_ - (X_max_ - X_min_) * 0.1

plt.plot([X_min , X_max],[0, 0], 'k-')

plt.ylim(y_min, y_max)
plt.xlim(X_min, X_max)
plt.title('distance from other sample vs std',fontsize=20)
plt.xlabel('distance from other sample',fontsize=20)
plt.ylabel('std',fontsize=20)

plt.show()
```
![2_distance_vs_std.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/348376/7c630512-603b-4ca1-eba4-cc18126054a5.png)

「distance from other sample」が 10以下の領域では平均距離が小さくなるにつれて、予測値のばらつきが小さくなっている（＝似たようなサンプルが多い領域での予測値は安定する）ことが確認できた。
10以上になるとばらつきはあまり平均距離に依存しなくなっている。



#参考HP
- [金子研究室HP](https://datachemeng.com/)
- [RDKitでケモインフォマティクスに入門(化学の新しいカタチ)](https://future-chem.com/rdkit-intro/)