予測モデルの適用範囲の可視化に関する一考察

はじめに

記事 論文と公共データベースを使って無料で始めるAI創薬　のおまけで予測モデルの適用範囲の可視化を行った。その際に、学習データのみで学習したPCA、UMAPによる次元圧縮モデルに対し、予測対象のデータを適用し可視化を行った。その際に一見、予測対象データが学習データの範囲内に入っているっぽいという結論を下したが、本当にそうなのか気になったので検証してみる。

検証内容

記事では、「図をみると予測対象データは、学習データの範囲から大きく逸脱していないように見えるが、学習データから構成されるデータのみで次元削減していることも影響しているかもしれない。」と書いた。
そこで、今回、学習データと予測対象データの全てを使って次元圧縮を行ってみて、学習データからのみ次元圧縮した場合の図と比較してみたい。

ソース

前回のソースに対し、学習データと予測対象データの全てを使って次元圧縮するように修正したソースを掲載する。

view_ad.py

import argparse
import csv

import pandas as pd
import numpy as np
import umap
import matplotlib.pyplot as plt
from sklearn.manifold import TSNE, MDS
from sklearn.decomposition import PCA

from rdkit import Chem
from rdkit.Chem import Descriptors, AllChem
from rdkit import rdBase, Chem, DataStructs

def main():

    parser = argparse.ArgumentParser()
    parser.add_argument("-train", type=str, required=True)
    parser.add_argument("-predict", type=str)
    parser.add_argument("-result", type=str)
    parser.add_argument("-method", type=str, default="PCA", choices=["PCA", "UMAP"])

    args = parser.parse_args()

    # all
    all_datas = []

    # all_train.csvの読み込み, fpの計算
    train_datas = []
    train_datas_active = []
    train_datas_inactive = []

    with open(args.train, "r") as f:
        reader = csv.DictReader(f)
        for row in reader:
            smiles = row["canonical_smiles"]

            mol = Chem.MolFromSmiles(smiles)
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=3, nBits=2048, useFeatures=False, useChirality=False)
            train_datas.append(fp)

            if int(row["outcome"]) == 1:
                 train_datas_active.append(fp)
            else:
                 train_datas_inactive.append(fp)

            all_datas.append(fp)



    if args.predict and args.result:
        result_outcomes = []
        result_ads = []

        # 予測結果読み込み
        with open(args.result, "r",encoding="utf-8_sig") as f:
            reader = csv.DictReader(f)
            for i, row in enumerate(reader):
                #print(row)
                if row["Prediction"] == "Active":
                    result_outcomes.append(1)
                else:
                    result_outcomes.append(0)

                result_ads.append(row["Confidence"])


        # drugbank.csvの読み込み, fpの計算
        predict_datas = []
        predict_datas_active = []
        predict_datas_inactive = []
        predict_ads = []
        with open(args.predict, "r") as f:
            reader = csv.DictReader(f)
            for i, row in enumerate(reader):
                print(i)
                smiles = row["smiles"]
                mol = Chem.MolFromSmiles(smiles)
                fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=3, nBits=2048, useFeatures=False, useChirality=False)
                predict_datas.append(fp)

                if result_outcomes[i] == 1:
                    predict_datas_active.append(fp)
                else:
                    predict_datas_inactive.append(fp)

                all_datas.append(fp)

    # 分析
    model = None
    if args.method == "PCA":
        model = PCA(n_components=2)
        model.fit(train_datas)
        #model.fit(all_datas)

    if args.method == "UMAP":
        model = umap.UMAP()
        model.fit(train_datas)
        #model.fit(all_datas)

    result_train = model.transform(train_datas)
    result_train_active = model.transform(train_datas_active)
    result_train_inactive = model.transform(train_datas_inactive)

    plt.title(args.method)
    #plt.scatter(result_train[:, 0], result_train[:, 1], c="blue", alpha=0.1, marker="o")
    plt.scatter(result_train_active[:, 0], result_train_active[:, 1], c="blue", alpha=0.5, marker="o")
    plt.scatter(result_train_inactive[:, 0], result_train_inactive[:, 1], c="blue", alpha=0.5, marker="x")

    # 予測(predict)
    if args.predict and args.result:

        result_predict = model.transform(predict_datas)
        result_predict_active = model.transform(predict_datas_active)
        result_predict_inactive = model.transform(predict_datas_inactive)

        #plt.scatter(result_predict[:, 0], result_predict[:, 1], c=result_ads, alpha=0.1, cmap='viridis_r')
        plt.scatter(result_predict_active[:, 0], result_predict_active[:, 1], c="red", alpha=0.1, marker="o")
        plt.scatter(result_predict_inactive[:, 0], result_predict_inactive[:, 1], c="red", alpha=0.1, marker="x")

    plt.show()


if __name__ == "__main__":
  import argparse
import csv

import pandas as pd
import numpy as np
import umap
import matplotlib.pyplot as plt
from sklearn.manifold import TSNE, MDS
from sklearn.decomposition import PCA

from rdkit import Chem
from rdkit.Chem import Descriptors, AllChem
from rdkit import rdBase, Chem, DataStructs

def main():

    parser = argparse.ArgumentParser()
    parser.add_argument("-train", type=str, required=True)
    parser.add_argument("-predict", type=str)
    parser.add_argument("-result", type=str)
    parser.add_argument("-method", type=str, default="PCA", choices=["PCA", "UMAP"])

    args = parser.parse_args()

    # all
    all_datas = []

    # all_train.csvの読み込み, fpの計算
    train_datas = []
    train_datas_active = []
    train_datas_inactive = []

    with open(args.train, "r") as f:
        reader = csv.DictReader(f)
        for row in reader:
            smiles = row["canonical_smiles"]

            mol = Chem.MolFromSmiles(smiles)
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=3, nBits=2048, useFeatures=False, useChirality=False)
            train_datas.append(fp)

            if int(row["outcome"]) == 1:
                 train_datas_active.append(fp)
            else:
                 train_datas_inactive.append(fp)

            all_datas.append(fp)



    if args.predict and args.result:
        result_outcomes = []
        result_ads = []

        # 予測結果読み込み
        with open(args.result, "r",encoding="utf-8_sig") as f:
            reader = csv.DictReader(f)
            for i, row in enumerate(reader):
                #print(row)
                if row["Prediction"] == "Active":
                    result_outcomes.append(1)
                else:
                    result_outcomes.append(0)

                result_ads.append(row["Confidence"])


        # drugbank.csvの読み込み, fpの計算
        predict_datas = []
        predict_datas_active = []
        predict_datas_inactive = []
        predict_ads = []
        with open(args.predict, "r") as f:
            reader = csv.DictReader(f)
            for i, row in enumerate(reader):
                print(i)
                smiles = row["smiles"]
                mol = Chem.MolFromSmiles(smiles)
                fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=3, nBits=2048, useFeatures=False, useChirality=False)
                predict_datas.append(fp)

                if result_outcomes[i] == 1:
                    predict_datas_active.append(fp)
                else:
                    predict_datas_inactive.append(fp)

                all_datas.append(fp)

    # 分析
    model = None
    if args.method == "PCA":
        model = PCA(n_components=2)
        #model.fit(train_datas)
        model.fit(all_datas)

    if args.method == "UMAP":
        model = umap.UMAP()
        #model.fit(train_datas)
        model.fit(all_datas)

    result_train = model.transform(train_datas)
    result_train_active = model.transform(train_datas_active)
    result_train_inactive = model.transform(train_datas_inactive)

    plt.title(args.method)
    #plt.scatter(result_train[:, 0], result_train[:, 1], c="blue", alpha=0.1, marker="o")
    plt.scatter(result_train_active[:, 0], result_train_active[:, 1], c="blue", alpha=0.5, marker="o")
    plt.scatter(result_train_inactive[:, 0], result_train_inactive[:, 1], c="blue", alpha=0.5, marker="x")

    # 予測(predict)
    if args.predict and args.result:

        result_predict = model.transform(predict_datas)
        result_predict_active = model.transform(predict_datas_active)
        result_predict_inactive = model.transform(predict_datas_inactive)

        #plt.scatter(result_predict[:, 0], result_predict[:, 1], c=result_ads, alpha=0.1, cmap='viridis_r')
        plt.scatter(result_predict_active[:, 0], result_predict_active[:, 1], c="red", alpha=0.1, marker="o")
        plt.scatter(result_predict_inactive[:, 0], result_predict_inactive[:, 1], c="red", alpha=0.1, marker="x")

    plt.show()


if __name__ == "__main__":
    main()

一応解説しておくと、プログラムの修正したところは、modelでfitする際の引数を、学習データから全データ（学習データ＋予測対象データ）に変更した点である。

     #model.fit(train_datas)
     model.fit(all_datas)

結果

一応、繰り返すと青が学習データで、赤が予測データとなる。重なって見にくい点はご容赦願いたい。

PCA

学習データのみでfittingした場合

学習データ＋予測対象データ全部でfittingした場合

UMAP

学習データのみでfittingした場合

学習データ＋予測対象データ全部でfittingした場合

考察

• PCA/UMAPともに、学習データのみでfittingした場合に、一見多くの予測対象データが適用範囲内にあるようにみえる。
• しかし学習データ＋予測対象データ全部でfittingした場合には、学習データの領域から大きく外側に存在するものが大量に存在する。
• つまり、前者の図を見て、予測対象データが適用領域にあるかどうかを判定するのは非常に危険であるといえる。
• なぜこのようなことが起こったかというと、前者の図は、あくまでも学習データ全体のみの傾向を考慮した次元圧縮モデルであり、その傾向に沿わないデータに対してfittingしても、モデルでうまくとらえられないためである。
• では、適用範囲の判断をどうしたらよいのか？ということであるが、1つには図に頼らずに学習データとの距離を測るための何等かの計算式を定義し、数値に頼るというとういことである。
• もう一つは、学習データと予測対象データの全体を含むより大きな化合物集合を用いて次元圧縮モデルを作成することである。
• 後者の場合、自然界に存在する全ての化合物を用意すればいいのではないか、と思うかもしれないが、その場合はデータ数が多くなりすぎて実質モデルが作れないこと、またあまりにも空間が大きすぎて、学習データと予測データを区別できない状況になってしまうという問題があるので、予測モデルのドメインに応じた適度な化合物集合を用意する必要があると考える（なんか論文ありそうだなー）。