読書メモ『Pythonで体感！医療とAIはじめの一歩』

概要

東京医科歯科大学が「データサイエンスリテラシー」の教育のために出版した書籍。
https://www.yodosha.co.jp/yodobook/book/9784758124188/

出版の意義は序文に明らか。

例えばX線写真について、X線に関する物理学的知識を全く持たず、ただそこに写った白い影が肺であるとして胸部X線写真を見ている医師と、なぜフィルム上に肺が白く写るのか理解したうえで読影している医師とで、そこから得られる情報に差が出ることは明らかでしょう。AIやデータサイエンスについても同じことが言えます。

なので、この本は「医療データを題材にした、機械学習の入門書」と捉えておくと良い。
後から内容を振りかえれるように軽くまとめ。

古典的機械学習・深層学習に共通で必要な概念

可視化

まずとにかくデータを可視化するところから始める。

前処理

データは適切に表現されているか、データに欠損が無いかなどを確認し、学習モデルに入力できるように整える工程。どちらかというと、実際のデータサイエンスではここが本番。

Encoding

カテゴリ変数を数値に符号化する工程。

多重共線性

重回帰において、説明変数どうしに強い相関がある場合、推定した係数が不安定になる現象。説明変数どうしの相関が無いように使う変数を限定するなど、前処理が必要。

<ここから自分調査>
数学的には、例えば線形重回帰 $y = {x^T}{\beta} + {\epsilon}$ ($x$と${\beta}$は特徴量次元のベクトル) では、最小二乗法で推定される${\beta}$は $y {\times} ({x}{x^T})^{-1}x$と計算されるが、この逆行列が計算できない(近似解法でも収束しづらくなる)ことによるもの......のはず。
https://ja.wikipedia.org/wiki/多重共線性
https://en.wikipedia.org/wiki/Ordinary_least_squares
<ここまで自分調査>

過学習の防止

与えられた学習データをすべて使うと、汎化性能がない(つまり、今後未知のデータが入力されたときに使い物にならない)モデルを作ってしまう可能性がある。
それを防止するために、例えばtrainデータとtestデータを分割して評価するホールドアウト法などが使われる。
特に、決定木・ランダムフォレストでは、設定を間違えるとtrain側のscoreが1.0になるので毎回気を付けて確認すること。

古典的機械学習

線形回帰

数値データの集合を、線形関数でフィッティングするアルゴリズム。

患者のBMIから、1年後の糖尿病の進行度スコアを予測する、線形単回帰モデルを作成。

import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import load_diabetes
raw = load_diabetes(as_frame=True, scaled=False)
#print(raw.DESCR)

x = raw.data[['bmi']]
y = raw.target

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2) #過学習防止

from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(x_train, y_train)
print(model.score(x_train, y_train)) #ここが1.0だったりtestより高すぎたり、過学習してないか要確認
print(model.score(x_test, y_test))

plt.figure()
plt.scatter(x, y, color='blue') #plotなら折れ線グラフ、scatterなら散布図
plt.plot(x, model.predict(x), color='red')
plt.title('Scatter: BMI vs. Diabetes possibility')
plt.xlabel('BMI')
plt.ylabel('Diabetes possibility')
plt.show()

ロジスティック回帰

正解が0か1の二値を取る場合に、シグモイド関数でフィッティングして二値分類(と各クラスにデータが属する確率の回帰)を行うアルゴリズム。

腫瘍の半径から、良性(1)か悪性(0)かを分類する、ロジスティック単回帰モデルを作成。

import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import load_breast_cancer
raw = load_breast_cancer(as_frame=True)
#print(raw.DESCR)

x = raw.data[['mean radius']]
y = raw.target #DSCRには212 - Malignant(悪性), 357 - Benign(良性)と書いてあるが、正規化されて良性=1、悪性=0になっている。

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2) #過学習防止

from sklearn.linear_model import LogisticRegression
model = LogisticRegression()
model.fit(x_train, y_train)
print(model.score(x_train, y_train))
print(model.score(x_test, y_test))

plt.figure()
plt.scatter(x, y, color='blue', label='raw data') #plotなら折れ線グラフ、scatterなら散布図
plt.scatter(x, model.predict_proba(x1d)[:, 1], color='red', label='predicted')  #ロジスティック回帰の場合、predictは[0,1]を返す。predict_probaを使うと、[クラス0である確率, クラス1である確率]を返す。
plt.title('Radius vs. Is-benign?')
plt.xlabel('Raduis')
plt.ylabel('Is-benign?')
plt.legend()
plt.show()

サポートベクトルマシン

サポートベクトル(決定境界に一番近いデータ)...と決定境界の距離を最大化することで、決定境界を学習する分類アルゴリズム。
データをより高次元に拡張して分類する「カーネルトリック」により、非線形データにも対応して分離性能を上げることが可能。
なお、SVMには特徴量の重要度の概念はない。重要度が必要な時には決定木・ランダムフォレストを使う。

腫瘍の10次元の特徴量をPCAで2次元に圧縮したうえで、その2次元の特徴量から、良性(1)か悪性(0)かを分類する、線形SVM分類モデルを作成。

<ここから自分調査>
PCAは、特徴量の次元圧縮に使われる手法の一つ。本書ではもともと10次元の特徴量をそのままSVMモデルに入力していたが、決定境界を可視化してみたかったため、PCAを介して次元圧縮するコードを追加してみた。
<ここまで自分調査>

import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import load_breast_cancer
raw = load_breast_cancer(as_frame=True)
#print(raw.DESCR)

x = raw.data.iloc[:, 0:10]   #パラメータは30個あるが、10個の特徴量がmeanとstdとworstで分けられたもの。多重共線性を回避するため、meanだけを使うことにする。
y = raw.target  #0 = 悪性、1 = 良性

#可視化したいため、まずはxを次元削減して2次元にする。ここではPCAを使用。
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
x_new = pca.fit_transform(x)

plt.figure()
plt.scatter(x_new[y == 0][:, 0], x_new[y == 0][:, 1], color='red', label='y=0, Malignant') # y=0の点（赤）
plt.scatter(x_new[y == 1][:, 0], x_new[y == 1][:, 1], color='blue', label='y=1, Benign') # y=1の点（青）
plt.title('PCA')
plt.xlabel('1st principal')
plt.ylabel('2nd principal')
plt.legend()
plt.show()

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x_new, y, test_size=0.2)  #過学習防止

from sklearn.svm import SVC #support vector classifier
model = SVC(kernel='linear')  #SVMでは、使うカーネル関数を選べる。linearとrbf(Radial basis function)を選べ、後者はカーネルトリックを使うことで非線形境界を学習可能。
model.fit(x_train, y_train)
print(model.score(x_train, y_train))
print(model.score(x_test, y_test))

from sklearn.inspection import DecisionBoundaryDisplay
disp = DecisionBoundaryDisplay.from_estimator(
    model, x_new, response_method="predict",
    xlabel='1st principal', ylabel='2nd principal', title='PCA', 
    alpha=0.5,
)
disp.ax_.scatter(x_new[y == 0][:, 0], x_new[y == 0][:, 1], color='red', label='y=0, Malignant') # y=0の点（赤）
disp.ax_.scatter(x_new[y == 1][:, 0], x_new[y == 1][:, 1], color='blue', label='y=1, Benign') # y=1の点（青）
plt.show()

決定木

データを条件によって選り分けていく樹形図を構築することで、クラス分類(や回帰)を行うアルゴリズム。可視化力が圧倒的。

腫瘍の10次元の特徴量から、良性(1)か悪性(0)かを分類する、決定木分類モデルを作成。

import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import load_breast_cancer
raw = load_breast_cancer(as_frame=True)
#print(raw)

x = raw.data.iloc[:, 0:10] #多重共線性を避けるためにmeanのみ使う
y = raw.target  #0: malignant 悪性 / 1: benign　良性

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2) #過学習防止

from sklearn.tree import DecisionTreeClassifier
model = DecisionTreeClassifier(
    criterion='gini', #純度評価法。'gini'と'entropy'を選べるが、大きな違いはないと本には書かれている。
    max_depth=3, #深さの最大値。これを指定しないと、すべてのデータが分類されるまで分割を繰り返すので、過学習する。
    min_samples_split=10,  #一度の分岐を作るのに必要な最低サンプル数。これが小さいと、やたら細かく分類し始めるので、過学習する。
)
model.fit(x_train, y_train)
print(model.score(x_train, y_train))
print(model.score(x_test, y_test))

from sklearn.tree import plot_tree
plt.figure(figsize=(20,10))
plot_tree(model,
          feature_names=raw.feature_names,  #特徴量の名前。これを指定すると、可視化した図の中に名前を含めてくれる。
          class_names=raw.target_names, #同じく、目的変数(分類クラス)の名前。
          filled=True)
plt.show()

決定木 × バギング型アンサンブル学習 = ランダムフォレスト

決定木は単体では学習精度が高くない「弱学習機」である。弱学習機を複数作成し、多数決を取ることで性能を上げ、かつ過学習を防ぐ「アンサンブル学習」の手法の一つに、「ランダムフォレスト」がある。

「バギング」とは、全てのデータから(重複を許しながら)一部のデータを取り出し、それぞれのデータで弱学習機を学習させるアンサンブル学習の手法のこと。このバギング型アンサンブル学習で、弱学習機に決定木を用いるのが、ランダムフォレストである。

腫瘍の10次元の特徴量から、良性(1)か悪性(0)かを分類する、ランダムフォレスト分類モデルを作成。

import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import load_breast_cancer
raw = load_breast_cancer(as_frame=True)
#print(raw.DESCR)

x = raw.data.iloc[:, 0:10] #多重共線性を避ける為にmeanのみ使う
y = raw.target #0:悪性, 1:良性

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2) #過学習防止

from sklearn.ensemble import RandomForestClassifier
model = RandomForestClassifier(
    n_estimators=100,  #構築する決定木の数。多すぎると過学習する。試行錯誤必要。
    max_depth=3,  #木の深さ。設定しないと過学習する。
    max_features=5  # 特徴量の数の上限。木ごとに使える特徴量の数を制限することで過学習を避ける。デフォルトは、特徴量の総数の平方根。今回10次元にしているので、少し多めにしている。
)
model.fit(x_train, y_train)
print(model.score(x_train, y_train))
print(model.score(x_test, y_test))

#ランダムフォレストは、決定木単体のような木の可視化はできないが、特徴量ごとの重要度を出せる。
#PFI (Permutation Feature Importance) という指標で、特徴量ごとに、それを消して学習した場合と消さないで学習した場合の予測精度の差分を、重要度として出力する。
from sklearn.inspection import permutation_importance
importance = permutation_importance(model, x_test, y_test, n_repeats=10)
forest_pfi = pd.Series(
    importance.importances_mean,
    index=raw.feature_names[0:10]
)
forest_pfi.plot.bar()
plt.show()

深層学習

複数のモデルを包含できる高次のニューラルネットワークを指定して、誤差逆伝播で層間の重みを変化させていくことで、ニューラルネットワーク自身が最適なモデルを設計していく。

腫瘍の10次元の特徴量から、良性(1)か悪性(0)かを分類する、3層全結合ニューラルネットワーク(パーセプトロン)モデルを作成。

import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import load_breast_cancer
raw = load_breast_cancer(as_frame=False) #numpy配列として読み込む
#print(raw.DESCR)

x = raw.data[:, 0:10]  #numpyなので、ilocではなく直接スライスする。 #多重共線性を避ける為にmeanのみ使う
y = raw.target #0:悪性, 1:良性

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2) #過学習防止

from keras.models import Sequential
from keras.layers import Dense
model = Sequential()
model.add(Dense(4, input_shape=(10, ), activation='relu'))  #Denseは全結合の意味。入出力の次元は今は適当に設定したが、本来はハイパーパラメータチューニングが必要。sklearn.model_selection.GridSearchCVなど
model.add(Dense(1, activation='sigmoid'))
model.compile(loss='binary_crossentropy', optimizer='Adam', metrics=['accuracy']) #学習の仕方の設計。2値分類なら2値交差エントロピー。
model.summary()
result = model.fit(x_train, y_train, batch_size=32, epochs=300)

plt.plot(result.history['accuracy'])
plt.xlabel('epochs [-]')
plt.ylabel('accuracy [-]')
plt.show()

_, ev_accuracy = model.evaluate(x_test, y_test)
print(ev_accuracy)

また、肺のX線画像(4096次元のベクトル)から肺炎(1)か健康(0)かを予測する、5層のニューラルネットワークを作成。

###前処理
#まず、フォルダ内のファイル名のリストを作成する。リスト内法表記
import os
list_healthy = [i 
                for i in os.listdir('/content/images/COVID-NORMAL/healthy') 
                if not i.startswith('.')]
list_covid19 = [i
                for i in os.listdir('/content/images/COVID-NORMAL/covid19')
                if not i.startswith('.')]
num_healthy = len(list_healthy)
num_covid19 = len(list_covid19)
num_all = num_healthy + num_covid19

#つぎに、ファイルを読み込んでベクトル化する。
#入力(x)について：画像の枚数はnum_all、画像サイズは64×64、1ピクセル当たりの次元数は1(グレースケール)とする。
#出力(y)について：画像の枚数はnum_all、出力は1(肺炎)、0(健康)でシグモイド関数で肺炎の確率を出力するものとする。
#このとき、xの次元は(num_all, 64, 64, 1), yの次元は(num_all, 1)となる。
import numpy as np
x = np.zeros((num_all, 64, 64, 1), dtype=float)
y = np.zeros((num_all, 1), dtype=int)

from keras.preprocessing.image import load_img
from keras.preprocessing.image import img_to_array

path_healthy = '/content/images/COVID-NORMAL/healthy'
for i in range(num_healthy):
  filename = f'{path_healthy}/{list_healthy[i]}'
  file = load_img(filename, color_mode='grayscale', target_size=(64, 64))
  x[i] = img_to_array(file) / 255  #正規化
  y[i] = 0 #健康

path_covid19 = '/content/images/COVID-NORMAL/covid19'
for i in range(num_covid19):
  filename = f'{path_covid19}/{list_covid19[i]}'
  file = load_img(filename, color_mode='grayscale', target_size=(64, 64))
  x[num_healthy + i] = img_to_array(file) / 255  #正規化
  y[num_healthy + i] = 1 #肺炎

#前処理最後に、ベクトルをシャッフルしておく
num_list = np.arange(num_all)
np.random.shuffle(num_list)
x = x[num_list]
y = y[num_list]


###モデル作成・学習
#深層学習のfitには、エポックごとに検証データを分割して過学習防止してくれるvalidation_split機能がある。
#そのため、train_test_splitで分けておく必要はない。
from keras.models import Sequential
from keras.layers import Dense, Flatten, Dropout
model = Sequential()
model.add(Flatten(input_shape=(64, 64, 1)))  #配列を1次元に変換する
model.add(Dense(512, activation='relu'))     #中間層
model.add(Dropout(0.5))                      #過学習防止
model.add(Dense(256, activation='relu'))     #中間層
model.add(Dropout(0.5))                      #過学習防止
model.add(Dense(128, activation='relu'))     #中間層
model.add(Dense(1, activation='sigmoid'))    #0か1を出力
model.compile(loss='binary_crossentropy', optimizer='Adam', metrics=['accuracy']) #学習の仕方の設計
model.summary()

result = model.fit(x, y, batch_size=32, epochs=300, validation_split=0.2) #validation_splitで過学習防止

import matplotlib.pyplot as plt
plt.plot(result.history['accuracy'], label='train_accuracy')
plt.plot(result.history['val_accuracy'], label='validation_accuracy')
plt.legend()
plt.xlabel('epochs [-]')
plt.ylabel('accuracy [-]')
plt.show()

test1 = img_to_array(load_img('/content/images/covid.jpg', color_mode='grayscale', target_size=(64, 64))) / 255
test2 = img_to_array(load_img('/content/images/NORMAL.jpg', color_mode='grayscale', target_size=(64, 64))) / 255
test = np.zeros((2, 64, 64, 1))
test[0] = test1
test[1] = test2
test_pred = model.predict(test)
print(test_pred)

AIと医療の未来の話

この本の重要な部分ですが、どう要約しようとしてもネタバレが強くなってしまうので、この読書メモでは記載しないことにしました。。。公共データベースの話や医療AIベンチャーの話など、今まで知らなかったことが多く載っていて、初学者の人にもとても参考になると思います。