機械学習初心者がMoAを予測してみた

はじめに

kaggleの過去のコンペ「Mechanisms of Action (MoA) Prediction」(https://www.kaggle.com/c/lish-moa)のデータセットで予測モデルを作成してみます。それなりに精度が出たら、kaggleにsubmitしてscoreを確認してみようと思います。
実行環境はMac (M1)とGoogle Colaboratoryです。

1.コンペの概要と評価方法

• コンペの概要

MoA (作用機序)とは薬が治療効果を及ぼす仕組みのことです。 投与した薬剤がどの標的分子(タンパク質)に作用してどんな影響を与えた結果、治療効果が得られるのかということを明らかにしたものです。
このコンペの目的は、薬剤で処理した100種類の細胞の生存率と、772種類の遺伝子の発現データを特徴量として206種類のMoAを予測するモデルを作成することです。良い予測モデルができれば、作用機序が不明な薬剤についても、この薬剤を投与後の細胞生存率の変化や遺伝子発現などのデータを利用してMoAが予測できるようになります。
では、どんな方法で予測するかですが、206種類の予測値を出力するため、ニューラルネットワークモデルを作ることにしました。

• モデルの評価方法とMoA予測値の提出方法

モデルの評価は下記の式で行います。

N: sig_id(薬剤の数)、M: MoA(targets)の数、ˆy: 予測データ、y: 正解データ

全ての予測値に対して正解との誤差を求め、平均を出しています。scoreが０に近いほど精度が良いということになります。この式で、作成したモデルの評価を行います。コードは下記の通りです。

def prob_minmax(score):
  return max([min([score,1-10**-15]),10**-15])

def calc_logloss(y_pred,y_gt):
  score=0
  for m in range(y_pred.shape[1]):
    for n in range(y_pred.shape[0]):
      score += y_gt[n,m]*np.log(prob_minmax(y_pred[n,m])) + (1-y_gt[n,m])*np.log(1-prob_minmax(y_pred[n,m]))
  score=-1*score/y_pred.shape[0]/y_pred.shape[1]
  return score

MoA予測値は、付与されたsubmission_sample様式に従って提出します。

2.データのダウンロードと確認

提供されたデータをダウンロードして内容を確認します。

# 必要なライブラリをインポート
import pandas as pd
import numpy as np
from keras.models import Sequential,Model
from keras.layers import Dense,Dropout,Input,concatenate
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from sklearn.metrics import log_loss
from iterstrat.ml_stratifiers import MultilabelStratifiedKFold
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import itertools

# データセットを読み込む
df_X_train = pd.read_csv('/**/MoA/train_features.csv')
df_X_test = pd.read_csv('/**/MoA/test_features.csv')
df_y_train = pd.read_csv('/**/MoA/train_targets_scored.csv')
df_y_train_ns = pd.read_csv('/**/MoA/train_targets_nonscored.csv')
df_ss = pd.read_csv('/**/MoA/sample_submission.csv')

• 説明変数(特徴量)ーtrain_features

カラムの情報は下記の通りです。

① sig_id: 実験した薬剤(名前を匿名化している)
② cp_type: trt_cpが実薬、ctl_vehicleが対照薬(コントロール)
③ cp_time: 薬効を測定した時間  24,48,72時間の3回
④ cp_dose: 薬剤の用量。D1とD2の2種類
⑤ g-: 遺伝子の発現量。772種類
⑥ c-: 薬剤投与後の細胞の生存率。100種類

train_features

• 目的変数ーtrain_targets, train_targets_nonscored

MoAのラベル。このコンペでは評価対象のtrain_targets(206ラベル)の他に、評価の対象にならないtrain_targets_nonscored(402ラベル)も付与されました。後者は予測しても評価対象にはならないので、後にスタッキング学習(目次4-3)に使います。

train_targets

3.機械学習用にデータの前処理

• 文字列の削除 : dropメソッドで学習に不要な文字列'sig_id'を削除します。

• カテゴリ変数に変換：特徴量に含まれる文字列'cp_type'と'cp_dose'は数値に変換する必要があります。ただし、trainデータとvalidデータで共通の文字列は同じ値になるように、まずカテゴリ変数に変換してから(下図参照)、get_dummies()関数で数値に変換します。

カテゴリ変数に変換するコード

• tensor型(TensorFrowで使える形)に変換：入力データ型をastype()メソッドでfloat32に変換します。

前処理をまとめたコードは下記の通りです。

X_train = df_X_train.drop(columns='sig_id')
X_test = df_X_test.drop(columns='sig_id')

categories_cp_type = set(X_train['cp_type'].unique().tolist() + X_test['cp_type'].unique().tolist())
categories_cp_dose = set(X_train['cp_dose'].unique().tolist() + X_test['cp_dose'].unique().tolist())
print(categories_cp_type)
print(categories_cp_dose)

X_train['cp_type'] = pd.Categorical(X_train['cp_type'], categories=categories_cp_type)
X_train['cp_dose'] = pd.Categorical(X_train['cp_dose'], categories=categories_cp_dose)
X_test['cp_type'] = pd.Categorical(X_test['cp_type'], categories=categories_cp_type)
X_test['cp_dose'] = pd.Categorical(X_test['cp_dose'], categories=categories_cp_dose)

print(X_train['cp_type'].dtypes)

# cp_type,cp_doseをOne-Hot_Encodingで変換
X_train = pd.get_dummies(X_train, columns=['cp_type','cp_dose']).astype('float32')
X_test = pd.get_dummies(X_test, columns=['cp_type','cp_dose']).astype('float32')

y_train = df_y_train.drop(columns='sig_id').astype('float32').values
y_train_ns = df_y_train_ns.drop(columns='sig_id').astype('float32').values

4.モデルの作成

1.baseline

まずはbaselineとなるモデルを作成します。
x_trainとy_trainを学習データと評価データに分けます。特徴量に偏りが出ないように、k-分割交差検証（K-fold Cross-validation)を使います。今回は５分割して5回学習させ、平均scoreを出します。

mskf = MultilabelStratifiedKFold(n_splits=5, random_state=0, shuffle=True)
train_score_list, valid_score_list, model_list=[],[],[]
for i, (train_idx, val_idx) in enumerate(mskf.split(X_train, y_train)):
    print(f'# no. {i} ####################################')
    X_train_k, y_train_k = X_train[train_idx, :], y_train[train_idx, :]
    X_val_k, y_val_k = X_train[val_idx, :], y_train[val_idx, :]

    model = Sequential()

    model.add(Dense(128, input_dim=X_train.shape[1], activation='relu'))
    model.add(Dropout(0.5))
    model.add(Dense(128, activation='relu'))
    model.add(Dense(206, activation='sigmoid'))  
    model.compile(loss='binary_crossentropy', optimizer='adam')

    # モデルを学習させる
    history = model.fit(X_train_k, y_train_k, epochs=50, batch_size=512, verbose=1, validation_data=(X_val_k, y_val_k) )

    # 予測値を出力する
    y_train_pred = model.predict(X_train_k)
    y_valid_pred = model.predict(X_val_k)

    train_score_list.append(calc_logloss(y_train_pred,y_train_k))
    valid_score_list.append(calc_logloss(y_valid_pred,y_val_k))
    model_list.append(model)

print(f'train score: {sum(train_score_list)/len(train_score_list)}, valid score: {sum(valid_score_list)/len(valid_score_list)}')

baselineモデルでは、train_scoreが 0.012501、varid_scoreが 0.015968でした。これをベースに、精度を上げる方法を考えます。

2.特徴量を作る

モデルの精度を上げるために特徴量を作ります。特徴量の作り方には'減らす'方法と、'増やす'方法があります。両方試してみます。

① 特徴量を減らすーPCA分析で次元削減
PCA分析で不要な情報を落とすとモデルがコンパクトになり、学習しやすくなります。今回の特徴量にはc(細胞の生存率)とg(遺伝子発現量)の2種類があるので、下記3パターンに分けて行います。

1. cの特徴量のみ

2. gの特徴量のみ

3. cとgをそれぞれ

まず、標準化を行います。次に'c-'(または'g-')の列のみ取り出してPCA分析を行います。そして元の特徴量からc列を削除しPCA分析した新しいc列を加えます。標準化とPCA分析はそれぞれscikit-learnのStandardScalerとPCAクラスを使います。下記はc列のみのコードです。(g列のコードはcをgに置き換えるだけなので省略します)

ss = StandardScaler()
ss.fit(X_train)
X_train　=　pd.DataFrame(ss.transform(X_train),columns=X_train.columns)
X_test　=　pd.DataFrame(ss.transform(X_test),columns=X_test.columns)

column_list　=　list(X_train.columns)
c_column_list　=　[col for col in column_list if 'c-' in col]
X_train_c = X_train[c_column_list]
X_test_c = X_test[c_column_list]

pca = PCA(n_components=60)
pca.fit(X_train_c)
X_train_c = pd.DataFrame(pca.transform(X_train_c),columns=['c-pca-'+str(col) for col in range(1,61)])
X_test_c = pd.DataFrame(pca.transform(X_test_c),columns=['c-pca-'+str(col) for col in range(1,61)])

X_train = X_train.drop(columns=c_column_list)
X_test = X_test.drop(columns=c_column_list)

X_train = pd.concat([X_train,X_train_c],axis=1)
X_test = pd.concat([X_test,X_test_c],axis=1)

② 特徴量を増やすー差分特徴量追加
特徴量同士の差をとり、新たな特徴量として加えます。ただし872個の特徴量全ての差分を取ると膨大な数になってしまうので、組み合わせの数を考慮して分散の大きい(=特徴量として重要である)方から1％分のみ使います。分散の大きい方から1％とは標準偏差値がいくつまでの特徴量を取れば良いのか？標準偏差のヒストグラムと累積比率のグラフを書いて確認してみます。

X_train.describe()

std=pd.DataFrame(X_train.describe().loc['std']['g-0':'c-99']).sort_values('std',ascending=False)
std

fig, ax1 = plt.subplots(figsize=(8, 6))
n, bins, patches = ax1.hist(std['std'],bins=30)
y2 = np.add.accumulate(n)/n.sum()
x2 = bins[1:]
ax2 = ax1.twinx()
ax2.plot(x2,y2,color='r')

ヒストグラムと累積比率

累積比率で0.99(上位から1％)のヒストグラムに含まれる特徴量(緑)を使えばよさそうです。次のコードで抽出します。

df_hist=pd.DataFrame({'n':n , 'y2':y2})
n_samples_cross=int(df_hist[df_hist['y2']>=0.99]['n'].sum())
features_cross=std.index[:n_samples_cross]
features_cross

特徴量を抽出できました。これら全ての組み合わせで差分を取ります。組み合わせを作るにはitertools.combinations()を使います。作った差分特徴量を元の特徴量に加えます。

X_train_ft1_ft2_diff_dict={}
for ft1,ft2 in itertools.combinations(features_cross,2):
  X_train_ft1_ft2=X_train[ft1]-X_train[ft2]
  X_train_ft1_ft2_diff_dict[f'{ft1}_{ft2}_diff']=X_train_ft1_ft2

X_train_ft1_ft2_diff=pd.DataFrame(X_train_ft1_ft2_diff_dict)

X_test_ft1_ft2_diff_dict={}
for ft1,ft2 in itertools.combinations(features_cross,2):
  X_test_ft1_ft2=X_test[ft1]-X_test[ft2]
  X_test_ft1_ft2_diff_dict[f'{ft1}_{ft2}_diff']=X_test_ft1_ft2

X_test_ft1_ft2_diff=pd.DataFrame(X_test_ft1_ft2_diff_dict)


X_train=pd.concat([X_train,X_train_ft1_ft2_diff],axis=1)
X_test=pd.concat([X_test,X_test_ft1_ft2_diff],axis=1)

3.nonscoredデータを用いてスタッキング

モデルの精度をさらに上げる方法として、アンサンブル学習が用いられます。アンサンブル学習とは個々に学習した複数のモデルを融合させる方法です。今回はニューラルネットワークを使っていること、nonscoredデータ(402ラベル)が付与されていることから、これらを生かしてスタッキングを試します。scored(206ラベル)を予測できるmodel 1と、nonscoredデータを予測できるように学習したmodel 2を融合して、さらに学習したmodel 3で最終的に予測値を出すことで、精度の向上が期待できます。
学習済みモデルを次のモデルに加える際には、予測値のみ出力させたいので、学習済みモデルの層の固定を行います。(下図) 

スタッキングのイメージ

model１の学習。5回の学習のうち、一番scoreが良かったモデルをスタッキングに使います。

mskf = MultilabelStratifiedKFold(n_splits=5, random_state=0, shuffle=True)
train_score_list, valid_score_list, model_list=[],[],[]
for i, (train_idx, val_idx) in enumerate(mskf.split(X_train, y_train)):
    print(f'# no. {i} ####################################')
    X_train_k, y_train_k = X_train[train_idx, :], y_train[train_idx, :]
    X_val_k, y_val_k = X_train[val_idx, :], y_train[val_idx, :]

    model = Sequential()

    model.add(Dense(128, input_dim=X_train.shape[1], activation='relu'))
    model.add(Dropout(0.5))
    model.add(Dense(128, activation='relu'))
    model.add(Dense(206, activation='sigmoid'))  
    model.compile(loss='binary_crossentropy', optimizer='adam')

    history = model.fit(X_train_k, y_train_k, epochs=50, batch_size=512, verbose=1, validation_data=(X_val_k, y_val_k) )

    y_train_pred = model.predict(X_train_k)
    y_valid_pred = model.predict(X_val_k)

    train_score_list.append(calc_logloss(y_train_pred,y_train_k))
    valid_score_list.append(calc_logloss(y_valid_pred,y_val_k))
    model_list.append(model)

print(f'train score: {sum(train_score_list)/len(train_score_list)}, valid score: {sum(valid_score_list)/len(valid_score_list)}')
print(f'best model no. is {np.argmin(valid_score_list)}')

# 5回の学習で一番良いスコアのモデルを取り出す
model1 = model_list[np.argmin(valid_score_list)]

model2の学習。model1と同様に、5回の学習のうち、一番scoreが良かったモデルをスタッキングに使います。

mskf = MultilabelStratifiedKFold(n_splits=5, random_state=0, shuffle=True)
train_score_list, valid_score_list, model_list=[],[],[]

for i, (train_idx, val_idx) in enumerate(mskf.split(X_train, y_train_ns)):
    print(f'# no. {i} ######################################')
    X_train_k, y_train_k = X_train[train_idx, :], y_train_ns[train_idx, :]
    X_val_k, y_val_k = X_train[val_idx, :], y_train_ns[val_idx, :]

    model = Sequential()

    model.add(Dense(128, input_dim=X_train.shape[1], activation='relu'))
    model.add(Dropout(0.5))
    model.add(Dense(128, activation='relu'))
    model.add(Dense(402, activation='sigmoid'))  
    model.compile(loss='binary_crossentropy', optimizer='adam')

    history = model.fit(X_train_k, y_train_k, epochs=35, batch_size=512, verbose=1, validation_data=(X_val_k, y_val_k) )

    y_train_pred = model.predict(X_train_k)
    y_valid_pred = model.predict(X_val_k)

    train_score_list.append(calc_logloss(y_train_pred,y_train_k))
    valid_score_list.append(calc_logloss(y_valid_pred,y_val_k))
    model_list.append(model)

print(f'train score: {sum(train_score_list)/len(train_score_list)}, valid score: {sum(valid_score_list)/len(valid_score_list)}')
print(f'best model no. is {np.argmin(valid_score_list)}')

# 5回の学習で一番良いスコアのモデルを取り出す
model2 = model_list[np.argmin(valid_score_list)]

model1とmodel2の層を固定し学習させないようにします。

for layer in model1.layers:
  layer.trainable = False

for layer in model2.layers:
  layer.trainable = False   

model3の学習。複数の入力層を作れるfunctional APIを使います。

mskf = MultilabelStratifiedKFold(n_splits=5, random_state=0, shuffle=True)
train_score_list, valid_score_list, model_list=[],[],[]
for i,(train_idx, val_idx) in enumerate(mskf.split(X_train, y_train)):
    print(f'# no.{i} #########################################')
    X_train_k, y_train_k = X_train[train_idx, :], y_train[train_idx, :]
    X_val_k, y_val_k = X_train[val_idx, :], y_train[val_idx, :]

    inputs = Input(shape=(None,X_train.shape[1]))
    x1 = model1(inputs)
    x2 = model2(inputs)
    x_concat = concatenate([x1,x2])
    x_concat = Dense(128,activation='relu')(x_concat)
    x_concat = Dense(128,activation='relu')(x_concat)
    x_concat = Dense(128,activation='relu')(x_concat)
    preds = Dense(206, activation='sigmoid')(x_concat)
    
    model_concat = Model(inputs=inputs, outputs=preds)
    model_concat.compile(loss='binary_crossentropy', optimizer='adam')

    history = model_concat.fit(X_train_k, y_train_k, epochs=70, batch_size=512, verbose=1, validation_data=(X_val_k, y_val_k) )

    y_train_pred = model_concat.predict(X_train_k)
    y_valid_pred = model_concat.predict(X_val_k)

    train_score_list.append(calc_logloss(y_train_pred,y_train_k))
    valid_score_list.append(calc_logloss(y_valid_pred,y_val_k))
    model_list.append(model_concat)

print(f'train score: {sum(train_score_list)/len(train_score_list)}, valid score: {sum(valid_score_list)/len(valid_score_list)}')
print(f'best model no. is{np.argmin(valid_score_list)}')
model3 = model_list[np.argmin(valid_score_list)]

5.scoreの比較とsubmit

ここまでで、モデルの精度を上げる準備ができました。PCA分析、差分特徴量追加、スタッキングのどの組み合わせで一番良いscoreを出せるでしょうか。実行した結果を表にしました。

c列のみPCA分析+差分特徴量を追加+スタッキングの組み合わせで最も良いscoreが出ました。

最後に、kaggleにsubmitしてみます。テスト用データ(X_test )を使ってこのモデルに予測させ、予測値をsubmitします。

y_test_pred = model3.predict(X_test)

df_ss.iloc[:,1:] = y_test_pred
df_ss.to_csv("submission.csv", index=False)

結果はPrivate Scoreが0.01730、Public Scoreが0.01966でした。だいぶ悪くなってしまいました。

全体を振り返ってみると、まだたくさん試せることがあります。まずパラメータチューニングを手動で行なっていたため、自動化できれば良かったかもしれません。また、アンサンブル学習においても、他の学習モデル(例えばc列とg列をそれぞれ学習させたものや、K-foldで5回に分けて学習したモデル全て)を使う、model3の学習時にさらに特徴量を加える、なども考えられます。

おわりに

現在、研究所で微生物の実験をしている私にとって、実験データから作用機序を予測するというこのテーマは大変興味深いものでした。今後、自分の業務に機械学習をどう取り入れられるか考えるモチベーションにもなりました。これからも機械学習の勉強を続けていきたいと思います。