交差検証を用いてElastic Netを化合物の溶解度データに対して最適化

「線形モデルを用いた化合物の溶解度予測：通常最小二乗法，Ridge回帰，Lasso回帰」という記事では，線形モデルと呼ばれる手法を用いて化合物の溶解度を予測する機械学習モデルを構築しました．

特に

特徴量が多い場合には通常最小二乗法では容易に訓練用データに過剰適合してしまうこと
過剰適合を抑えるための正則化法にはRidge回帰とLasso回帰という方法があること
正則化の程度を制御するにはαというハイパーパラメーターを変化させること

などを学びました．今回はRidgeとLassoの正則化の方法を組み合わせたElastic Netという手法を用いて，より効率的で信頼性の高いモデルの構築方法について学んでいきたいと思います．

今回も下記論文に付属のデータセットを用いて機械学習を行っていきます．
「ESOL: Estimating Aqueous Solubility Directly from Molecular Structure」
J. Chem. Inf. Comput. Sci. 2004, 44, 1000–1005.
(DOI: 10.1021/ci034243x)

Elastic Netとは
交差検証とは
1. テスト用データからの情報のリーク
グリッドサーチとは
scikit-learnを用いた交差検証とグリッドサーチ
終わりに

Elastic Netとは

Elasttic NetはRidge回帰によるL2正則化と，Lasso回帰によるL1正則化を合わせたモデルになります．コスト関数としては

正則化の程度を示すα
L1正則化とL2正則化の混合比を示すr

の2つをハイパーパラメーターとして，下のように示すことができます．

$$ Cost = MSE(w) + r \alpha \sum_{i=1}^{n} |w_{i}| + \frac{1-r}{2} \alpha \sum_{i=1}^{n} w_{i}^2 $$

Elastic Netとはscikit-learnではsklearn.linear_modelに実装されています．

交差検証とは

テスト用データからの情報のリーク

例えば「線形モデルを用いた化合物の溶解度予測：通常最小二乗法，Ridge回帰，Lasso回帰」という記事では，Ridge回帰の正則化の程度を制御するパラメータαの値を

あるαについて訓練用データを用いてモデル構築
作成したモデルのテスト用データに対する精度をチェック
得られた精度をもとに最適なαの値を決定

といった順番で評価しました．コードの一部を以下に抜粋します（注：このままでは動きません）．

for i in [0.001, 0.1, 1, 10, 100]:
    r = Ridge(alpha=i).fit(X_train, y_train)
    print('-- alpha = {} --'.format(i))
    print('acc on train: {:.3f}'.format(r.score(X_train, y_train)))
    print('acc on test: {:.3f}'.format(r.score(X_test, y_test)))

-- alpha = 0.001 --
acc on train: 0.983
acc on test: 0.231
-- alpha = 0.1 --
acc on train: 0.981
acc on test: 0.493
-- alpha = 1 --
acc on train: 0.959
acc on test: 0.637
-- alpha = 10 --
acc on train: 0.859
acc on test: 0.695
-- alpha = 100 --
acc on train: 0.613
acc on test: 0.570

得られた結果から，α=10の時にテスト用データに対するスコアが0.695と最も高いです．それではこの結果をもって，未知のデータに対しても同様の精度が期待できると言えるでしょうか？

残念ながら未知データに対する精度はテスト用データに対するものよりも低くなります．その理由としては，今回のように

訓練用データでモデル構築を行い
テスト用データの精度を毎回ハイパーパラメーターを変えるごとにチェック

という手順を踏むことで，テスト用データがハイパーパラメーターの最適化に使われてしまったためです．これをテスト用データの漏洩，リークといいます．

リークを防ぐためにはハイパーパラメーターをテスト用データを使わずに決める必要があります．そのためには，訓練用データをさらに「訓練用データ」と「評価用データ」に分けてパラメータの最適化を行うことが必要です．例えばテスト用データを除いた上で，残りの8割を訓練用に，2割を評価用に分けるとすると下図のようになります．

このように訓練用データと評価用データに分けてハイパーパラメーターを決定することでテスト用データの漏洩を防ぐことができます．欠点としては上図からも明らかなように，モデルの訓練に使えるデータ数が減ってしまう点が挙げられます．一般的にデータ数を多くしてモデルを構築した方が精度がよくなりますので，これは問題です．

またランダムに訓練用データと評価用データを分けることから，「たまたま訓練用データに学習の容易なデータが入り込んだ」場合には，評価用データに対する精度は低くなるでしょう．逆の可能性もありえます．

これらの問題に対する解決策として行われるのが「交差検証」（クロスバリデーション）と呼ばれる方法になります．交差検証を用いると単に１回だけ訓練用と評価用データに分割するよりも信頼性の高いモデルが得られます．

よく使われるのがk-分割交差検証（k-fold crossvalidation）と呼ばれるものです．これはテスト用データを除いた残りのデータをk個に分割し，1つのデータ群を評価用に，残りを訓練用に用います．これを評価用データの集団を変えながらk回繰り返します．下には5個に分割した場合を示しています．

5-foldの交差検証の結果，1つのモデル（パラメータ）に対して5個の精度が手に入ります．その平均値がモデルの訓練用データにおける精度と言えます．またばらつきを見ることで，構築したモデルの剛健性を評価することもできます．これは1回だけ訓練する場合と異なる点です．

また交差検証では結果として全てのデータが1回評価用データに入りますので，上で述べたような「データ分割におけるたまたま」の要素を排除でき，全てのデータに対してきちんと学習できないと評価が上がりません．さらに分割数を大きくすることで訓練用データの割合を増やせますので，効率的なデータの利用が可能です．

k分割交差検証を極限にまで突き詰めると，1つのサンプルを検証用とし残りのデータを訓練用とする分割法に行き着きます．この方法は「1つ抜き交差検証」（leave-one-out）と呼ばれています．1つ抜き交差検証は特にサンプル数が少ない場合に有効です．

当然ですが，k分割交差検証でkを大きくしたり，1つ抜き交差検証を用いることで構築するモデルの数が増えますので，より長い計算時間が必要となる点が欠点になります．自分が利用できる計算資源と時間を考えながら，モデル構築を行う必要があります．

pythonの機会学習用ライブラリであるscikit-learnでは「k分割交差検証」「1つ抜き交差検証」ともにsklearn.model_selectionに実装されています．具体的な使い方については後ほど見ていきます．

グリッドサーチとは

複数のハイパーパラメーターをグリッド上に変化させながら，最適な値を探索していく手法をグリッドサーチと呼びます．

今回扱うElastic Netの場合，

正則化の程度を示すα
L1正則化とL2正則化の混合比を示すr

の2つがハイパーパラメーターですので，2次元のグリッドで探索していきます．例えばαをm個，rをn個の値試す場合，下記のようにm・n個の精度が得られます．

	α₁	α₂	・・・	α_m
r₁	x_1,1	x_2,1	・・・	x_m,1
r₂	x_1,2	x_2,2	・・・	x_m,2
・・・
r_n	x_1,n	x_2,n	・・・	x_m,n

各パラメータのペアに対してk分割交差検証を行う場合，モデル構築はk・m・n回行われることになります．

scikit-learnを用いた交差検証とグリッドサーチ

それでは具体的にscikit-learnを用いた交差検証の使い方を見ていきましょう．最初にRDKitを用いて必要な分子の準備を行います．

必要なライブラリのインポートと分子の読み込み

まずはデータセットから分子を読み込み，分子の情報を入力ベクトルとして変換しましょう．

論文のデータはコンマ区切りのテキストファイルで，「化合物名」「溶解度の実験値」「（論文中の）溶解度の予測値」「SMILES」の4つのデータが1144化合物について格納されています．

以下のコードでは

必要なライブラリのインポート
pandasのデータフレームとして情報を読み込む
SMILESからRDKitのMOLオブジェクトを構築
Morganフィンガープリントを入力ベクトルとして準備
データセットを訓練用とテスト用に分割

という処理を行っています．

RDKitとpandasの併用方法については，「RDKitのPandasToolsでデータ分析を加速する」という記事で解説しています．参照してみてください．

## 1. ライブラリのインポート
from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, PandasTools
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib as mpl
import matplotlib.pyplot as plt
print(rdBase.rdkitVersion) # 2018.09.2
mpl.style.use('seaborn')
 
## 2. データフレームとして読み込み
df = pd.read_csv('./ci034243xsi20040112_053635.txt')
## 3. MOLオブジェクトの構築
PandasTools.AddMoleculeColumnToFrame(frame=df, smilesCol='SMILES')
 
## 4. フィンガープリントの作成
fps = []
for mol in df.ROMol:
    fp = [x for x in AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048)]
    fps.append(fp)
fps = np.array(fps)
morgan_fps.shape # (1144, 2048)
targes = df['measured log(solubility:mol/L)']
targes.shape # (1144,)
 
## 5. データセットの分割
from sklearn.model_selection import train_test_split
X_trainval, X_test, y_trainval, y_test = train_test_split(fps, targets, random_state=0)
print(X_trainval.shape, X_test.shape, y_trainval.shape, y_test.shape)
# ((858, 2048), (286, 2048), (858,), (286,))

scikit-learnを用いたk分割交差検証

linear_model.ElasticNet(alpha, l1_ratio)
model_selection.cross_val_score(model, data, cv)

k分割交差検証はcross_val_scoreで行うことができます．パラメータcvに数字を与えることで分割数を指定します．下の例では試しにα=0.01, r=0.5のElastic Netでモデル構築を行っています．l1_ratioがrに相当します．

from sklearn.linear_model import ElasticNet
from sklearn.model_selection import cross_val_score
elastic = ElasticNet(alpha=0.01, l1_ratio=0.5, max_iter=100000)
scores = cross_val_score(elastic, X_trainval, y_trainval, cv=5, n_jobs=-1)
print(scores)
print(scores.mean())

得られた結果はnumpy arrayとして格納されています．5個のスコアが得られています．

[0.64337466 0.61231278 0.63217715 0.71466749 0.62917055]
0.6463405239222204

scikit-learnを用いたグリッドサーチ

model_selection.GridSearchCV(model, parm_grid, cv)

scikit-learnではGridSearchCVを用いることでグリッドサーチを容易に行うことが可能です．探索するパラメータの組み合わせはparam_gridで指定します．

下記のコードではalphaを6種類，l1_ratioを5種類指定しています．各組み合わせに10-分割交差検証を用いていますので全部で300回のモデル構築を行います．今回はn_jobs=-1を指定することで利用可能なコアを全て用いています．作成したGridSearchCVオブジェクトはscikit-learnの他のモデルと同じようにfitメソッドで訓練を行うことができます．

なお訓練したモデルの中で最良の結果を与えたものはbest_estimator_属性でアクセス可能です．

from sklearn.model_selection import GridSearchCV
param_grid = {'alpha': [0.001, 0.01, 0.1, 1, 10, 100],
              'l1_ratio': [0, 0.25, 0.5, 0.75, 1]}
grid_search = GridSearchCV(ElasticNet(max_iter=100000), param_grid, cv=10, n_jobs=-1)
grid_search.fit(X_trainval, y_trainval)

grid_search.best_estimator_

ElasticNet(alpha=0.01, copy_X=True, fit_intercept=True, l1_ratio=0,
      max_iter=100000, normalize=False, positive=False, precompute=False,
      random_state=None, selection='cyclic', tol=0.0001, warm_start=False)

またグリッドサーチにおける各結果はcv_results_属性で取得できます．そのままでは一覧性が低いため，pandasのデータフレームとして保存しています．

下のコードではさらにデータフレームから交差検証におけるスコアの平均値を取得し，alphaとl1_ratioの値ごとに並べた6 x 5のnumpy配列を得ています．

grid_scores = pd.DataFrame(grid_search.cv_results_)
scores = np.array(grid_scores.mean_test_score).reshape(6,5)
print(scores)

[[ 0.60755904  0.61233831  0.61187027  0.61753221  0.6209419 ]
 [ 0.67229286  0.66291973  0.64676282  0.62781539  0.61138877]
 [ 0.54591023  0.41330003  0.31931601  0.24652866  0.19535217]
 [ 0.2268953   0.01075187 -0.02161896 -0.02161896 -0.02161896]
 [ 0.02278557 -0.02161896 -0.02161896 -0.02161896 -0.02161896]
 [-0.01673102 -0.02161896 -0.02161896 -0.02161896 -0.02161896]]

グリッド毎に可視化してみると，

l1_ratioによらず大きいalpha値ではスコアが低いこと
best_estimator_で見たようにalpha=0.01, l1_ratio=0が最良であること

が確認できます．

fig = plt.figure()
ax = fig.add_subplot(111)
mat = ax.matshow(s, cmap='Reds', alpha=0.8)
fig.colorbar(mat)
ax.set_xticks(range(5))
ax.set_xticklabels([0, 0.25, 0.5, 0.75, 1])
ax.set_xlabel('alpha')
ax.set_yticks(range(6))
ax.set_yticklabels([0.001, 0.01, 0.1, 1, 10, 100])
ax.set_ylabel('l1 ratio')

scikit-learnのElasticNetCVクラスを用いたグリッドサーチ

linear_model.ElasticNetCV(alphas, l1_ratio, cv)

上述したGridSearchCVを用いたハイパーパラメーターのスクリーニングはどのような機械学習モデルに対しても利用可能な一般性のある方法です．

一方でscikit-learnには一部の機械学習アルゴリズムに対しては自動で交差検証を用いたハイパーパラメーターのスクリーニングを行うクラスが実装されており，ElasticNetにもそういったものが存在します．

ElasticNetCVにスクリーニングしたいalphとl1_ratioをリストで与えることで交差検証を行い，最適なハイパーパラメーターを選びます．

下のコードでは既にalpha=0.01, l1_ratio=0近辺が最適とわかっていますので，この周辺をさらに細かくスクリーニングしています．alphaとして5個，l1_ratioを2個設定し，cv=10で交差検証を行います．

得られたベストのパラメータはそれぞれalpha_とl1_ratio_でアクセス可能です．また自動的にこれらのパラメータを用いて全てのデータを用いて訓練したモデルが作成されます．

from sklearn.linear_model import ElasticNetCV
elastic_cv = ElasticNetCV(alphas=[0.005, 0.01, 0.025, 0.05, 0.075], 
                          l1_ratio=[0, 0.025], 
                          max_iter=100000, cv=10, n_jobs=-1)
elastic_cv.fit(X_trainval, y_trainval)
print(elastic_cv.alpha_, elastic_cv.l1_ratio_)
print('trainval score: {:.3f}'.format(elastic_cv.score(X_trainval, y_trainval)))
print('test score: {:.3f}'.format(elastic_cv.score(X_test, y_test)))

0.075 0.0
trainval score: 0.674
test score: 0.613

今回の場合l1_ratio=0ですから，実質的にL2正則化のRidge回帰モデルを作成していることになります．「線形モデルを用いた化合物の溶解度予測：通常最小二乗法，Ridge回帰，Lasso回帰」という記事でも交差検証を用いずにRidge回帰によるalphaの値をスクリーニングし，alpha=10が最もよいスコアを与えました．

冒頭でも述べたように，以前の方法ではテスト用データの漏洩が起きていますので，未知データに対する予測精度としては今回のモデルの方が信頼のおける結果と言えます．

終わりに

今回は「交差検証を用いてElastic Netを化合物の溶解度データに対して最適化」という話題について，

テスト用データの漏洩と交差検証とはなにか
グリッドサーチによるハイパーパラメーターの最適化
Elastic Netを用いた交差検証の具体例

などを学んできました．特にテスト用データのリークは非常に重要な内容ですので，きちんと理解しておきましょう．

これまで「RDKitとscikit-learnで機械学習：変異原性をk-最近傍法で予測」という記事から数回にわたり，機械学習の基本から，いくつかのアルゴリズムを用いて化合物の変異原性の有無や溶解度の予測を行うモデルを構築してきました．その際，分子の入力情報としては構造から得られるフィンガープリントを用いてきました．すなわち0か1の値から構成される疎なベクトルとして分子を表現してきました．

本サイトでもさまざまな分子記述子を扱ってきたように，入力情報の表現方法としては他にも色々考えられます．たとえばリピンスキールールのように「分子量」「脂溶性（CLogP）」「水素結合受容体・供与体の数」を基にして入力ベクトルを作成することも可能でしょう．

次回はこういった分子記述子を用いて分子の入力ベクトルを作成する場合に注意したほうがよいデータの前処理について，「k-最近傍法」や「サポートベクターマシン」というアルゴリズムを例にして見ていきたいと思います．

＞＞次の記事：「データ分析と前処理：パイプライン処理で化合物の溶解度を推定」