背景
RDKitやmordred等記述子計算で生成した特徴(記述計算結果)を元に、化合物の機械学習を行う際、そのままでは特徴の数が多くオーバーフィッティングの可能性があるため、特徴選択を行う必要がある。今回、基本的な特徴選択アルゴリズムであるscikit-learnのVarianceThesholdを試したので記載する。
環境
- Windows10
- RDkit 2018.09.1.0
- scikit-learn 0.20.2
- Jupyter Notebook
VarianceThresholdとは
参考文献[1]を読んだ内容を簡単にまとめると、所定の分散を満たさないサンプルを持つ特徴を除去するものである。デフォルトでは0分散つまり全てのサンプルが同じ値である特徴が除去される。他の例として, 0, 1のブール値を持つ特徴に対し、サンプルの80%以上が0または1の特徴を削除する例が記載されている。この場合サンプルが、0,1の値がベルヌーイ確率変数に従うと仮定した場合、分散は0.8*(1-0.8)であるから、これを引数に与えればよい。
確かにどのデータもほとんど同じであるような特徴を使っても、目的変数の説明には役立ちそうにないから、除外することは理にかなっているといえる。
試したこと
以下Jupyter NoteBookで試したことを順に説明する。
1. 必要データのimport
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
from rdkit import Chem
from rdkit.Chem import Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors
%matplotlib inline
2. データの読み込み
今回、MoleculeNetから取得できるBACEのデータ(1513化合物)を使った。このデータは
https://s3-us-west-1.amazonaws.com/deepchem.io/datasets/molnet_publish/bace.zip
からダウンロードできる。解凍しcsvをパスを指定してpandasのDataFrame読み込む。
df = pd.read_csv("../data/bace/bace.csv")
3. SMILESをRDKitのMOLオブジェクトに変換
pandasのDataFrameを見ると、0番目がsmilesとなっているため、このデータを取り出し、
RDKitでMOLオブジェクトに変換する。
smiles_all = df.iloc[:, 0].values
mols = []
for smi in smiles_all:
mol = Chem.MolFromSmiles(smi)
if mol != None:
mols.append(mol)
else:
print(smi)
4. RDKitで記述子計算
RDKitで記述子計算を行う。
descs = [desc_name[0] for desc_name in Descriptors._descList]
desc_calc = MoleculeDescriptors.MolecularDescriptorCalculator(descs)
result = []
for i, mol in enumerate(mols):
result.append(desc_calc.CalcDescriptors(mol))
desc_df = pd.DataFrame(result)
desc_df.columns = descs
desc_df.head()
desc_dfに200個の特徴が生成されている。
5. VarianceThresholdで特徴選択
以下のようにデフォルトの引数でVarianceThresholdを実行し、
特徴選択を行った。
from sklearn.feature_selection import VarianceThreshold
X = desc_df.values
select = VarianceThreshold()
X_new = select.fit_transform(X)
その結果200個の特徴が168個に削減された。
試しに、どの特徴が除去されたか確認する。
np.array(descs)[select.get_support()==False]
array(['NumRadicalElectrons', 'SMR_VSA8', 'SlogP_VSA9', 'VSA_EState1',
'VSA_EState2', 'VSA_EState3', 'VSA_EState4', 'VSA_EState5',
'VSA_EState6', 'VSA_EState7', 'fr_C_S', 'fr_HOCCN', 'fr_SH',
'fr_aldehyde', 'fr_azide', 'fr_azo', 'fr_barbitur',
'fr_benzodiazepine', 'fr_diazo', 'fr_dihydropyridine',
'fr_epoxide', 'fr_hdrzine', 'fr_hdrzone', 'fr_imide', 'fr_isocyan',
'fr_isothiocyan', 'fr_lactam', 'fr_nitroso', 'fr_phos_acid',
'fr_phos_ester', 'fr_prisulfonamd', 'fr_thiocyan'], dtype='<U24')
例えば、上記にある"SMR_VSA8'のヒストグラムを表示してみる。
plt.hist(desc_df['SMR_VSA8'].values)
ちなみにsckit-learnでの特徴選択の結果は、get_support()メソッドにTrue, Falseのリストの形で格納されているため、それを使って、以下のように新たな特徴セットを作成することができる。
desc_new_df = pd.DataFrame(X_new, columns=np.array(descs)[select.get_support()==True])