【初心者向け】基礎&実践プログラミング

初心者がつまづきやすいところ、最短で実力が身につく方法をお伝えします。

【機械学習】ランダムフォレスト (Random Forest)で疾患分類

f:id:AIProgrammer:20200428214520p:plain (Credit: Decision tree icon Royalty Free Vector Image - VectorStock )

目的

自分のメモ用。 ランダムフォレスト (Random Forest)を使って疾患Aと疾患Bを分ける。
プログラミング言語はPythonで、ライブラリはscikit-learnを用いる。

入力したバイオマーカの内どれが重要かを調べる。 UCLの先行研究では、Accuracy74.0%だったのでそれを超えることが目標。

作業環境

$ cat /etc/os-release
NAME="Ubuntu"
VERSION="18.04.3 LTS (Bionic Beaver)"
ID=ubuntu
ID_LIKE=debian
PRETTY_NAME="Ubuntu 18.04.3 LTS"
VERSION_ID="18.04"
HOME_URL="https://www.ubuntu.com/"
SUPPORT_URL="https://help.ubuntu.com/"
BUG_REPORT_URL="https://bugs.launchpad.net/ubuntu/"
PRIVACY_POLICY_URL="https://www.ubuntu.com/legal/terms-and-policies/privacy-policy"
VERSION_CODENAME=bionic
UBUNTU_CODENAME=bionic

実行

主要なmoduleのimport

import modules
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib import pyplot as plt
import os
import numpy as np

Future warningがうざいのでここで無視するよう指示。

# ignore warning
import warnings
warnings.simplefilter('ignore')

resultをRF_resultに保存するため、RF_resultフォルダを作成。

# make directory
if not os.path.isdir('RF_result'):
    os.mkdir('RF_result')

データセットの読み込み。sys.argvを使うことで引数を指定してあげることができる。 ここでの引数はデータセットであるCSVファイルを指定してやればいい。 このデータセットにはコントロール群と疾患A, Bのデータがある。

print('Start Random Forest...')
# import csv file
import sys
args = sys.argv
import pandas as pd
df = pd.read_csv(args[1])

.str.contains()で利用したいバイオマーカを絞る。 そしてデータだけ取り出し。

# select biomarker
#variable = df.iloc[:,0:672] # variable is from 1st to 672nd column.
## select map
variable = df.loc[:,df.columns.str.contains('ICVF|OD|ISO|MVF')]

欠損値medianでfillする。

# processing for defect values.
variable = variable.fillna(variable.median()) 

#カテゴリ変数の変換
#df['Sex'] = df['Sex'].apply(lambda x: 1 if x == 'male' else 0)
#df['Embarked'] = df['Embarked'].map( {'S': 0, 'C': 1, 'Q': 2} ).astype(int)
#df = df.drop(['Cabin','Name','PassengerId','Ticket'],axis=1)

コントロール群であってもAge,Sex,Intracranial Volume(ICV)によって個人差がある。 これらの影響を除くための処理。 線形近似し、y=ax+bの傾きaをゼロにするようにデータを補正。

# adjust biomarkers for covariates
from sklearn import linear_model
def nu_freq(X_reg, y_reg):
    lr = linear_model.LinearRegression()
    lr.fit(X_reg, y_reg)
    return lr.coef_, lr.intercept_

def control_residuals(X, y, covar):
    control_covar = covar[y == 1].astype(float)
    control_X = X[y == 1].astype(float)
    n_biomarkers = control_X.shape[1]
    regr_params = dict((x, None) for x in range(n_biomarkers))
    for i in range(n_biomarkers):
        mask = ~np.isnan(control_X[:, i]) & ~np.isnan(control_covar).any(axis=1)
        X_reg, y_reg = control_covar[mask, :], control_X[mask, i]
        regr_params[i] = nu_freq(X_reg, y_reg)
    residuals = X.copy()
    for key, value in regr_params.items():
        residuals[:, key] -= regr_params[key][1]
        residuals[:, key] -= np.matmul(covar, regr_params[key][0])
    return residuals
## read in covariates
covars = df[['Age', 'Gender', 'ICV']].astype(float).values
## adjust biomarkers for covariates
adjusted_variable = control_residuals(variable.values[:,0:672], df.Groups.values, covars)
## save data 
variable.to_csv("./RF_result/pre_adjusted_variable.csv", index=False)
adjusted_variable_save = pd.concat([pd.DataFrame(variable.columns).T, pd.DataFrame(adjusted_variable)], axis=0)
adjusted_variable_save.to_csv("./RF_result/post_adjusted_variable.csv", index=False)

実際にやりたいのは、コントロール群、疾患A, Bの3分類ではなく、疾患A, Bの2分類なので 疾患A, Bのみをpick up。

## select subject
adjusted_variable = adjusted_variable_save[(df['Groups'] == 2)| (df['Groups'] == 3)]
#print(adjusted_variable)
Groups = df.Groups[(df['Groups'] == 2)| (df['Groups'] == 3)]

ついにRandom Forest(ランダムフォレスト)を実行。 ハイパーパラメータの最適化もしたいのでGrid search(グリッドサーチ)を実行。 (都合上最適化されたものを記載。)

# random forest
#from sklearn.ensemble import RandomForestClassifier
#from sklearn.metrics import (roc_curve, auc, accuracy_score)
#clf = RandomForestClassifier(n_estimators=100, random_state=0)
#clf = clf.fit(train_X, train_y) # training
#pred_y = clf.predict(test_X) # validation
#fpr, tpr, thresholds = roc_curve(test_y, pred_y, pos_label=1)
#print('AUC:', auc(fpr, tpr))
#print('Accuracy:', accuracy_score(pred_y, test_y))

# Grid Search
# Create the parameter grid based on the results of random search
## Create a based model
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV
param_grid = {
    'bootstrap': [True],
    'max_depth': [120], # [110,115,120,125,130]
    'max_features': [200], # [200,212,225,273,250]
    'min_samples_leaf': [4], # [4, 5]
    'min_samples_split': [8], # [8, 10, 12]
    'n_estimators': [130] # [100,110,120,130]
} 
clf = RandomForestClassifier(random_state=5) # random_state=5, 10, Accuracy:0.754 # Instantiate the grid search model
grid_search = GridSearchCV(estimator = clf, param_grid = param_grid, 
                          cv = 5, n_jobs = -1, verbose = 1)
## Fit the grid search to the data
grid_search = grid_search.fit(adjusted_variable, Groups) # training
print("\nBest Parameter is below.")
print(grid_search.best_params_)
print("\n")

Grid searchで最適化されたパラメータでone leave out(one-leave-out)を実行。 1件1件の正誤と平均のAccuracyが算出される。 今回の場合だと、75.4%。

# validation
from sklearn.model_selection import cross_val_score
from sklearn.metrics import recall_score
best_grid = grid_search.best_estimator_
## k-fold cross validation
#scores = cross_val_score(best_grid, adjusted_variable, Groups, cv=5)
## leave-one-out
from sklearn.model_selection import LeaveOneOut
scores = cross_val_score(best_grid, adjusted_variable, Groups, cv = LeaveOneOut()) #分類器としてLeaveOneOut()を指定
print('Cross-Validation scores: {}'.format(scores))
print('Average score: {}'.format(np.mean(scores)))

重要度(importance)の算出と図示。 さらに保存。

# importance
## visualyze
features = variable.columns
importances = best_grid.feature_importances_
indices = np.argsort(importances)
plt.figure(figsize=(10,10))
plt.barh(range(len(indices)), importances[indices], color='b', align='center')
plt.yticks(range(len(indices)), features[indices])
plt.savefig('RF_result/importance.png')

## save importance values
importance = pd.DataFrame({ 'Valiable' :features[indices], 'Importance' :importances[indices]})
importance.to_csv("./RF_result/importance.csv", index=False)

使用した決定木の保存。

# save the trees(forest)
import pydotplus
from IPython.display import Image
from graphviz import Digraph
from sklearn.externals.six import StringIO
from sklearn import tree
print('\nSaving result to "RF_result" folder...')
for i in range(1, 101):
    dot_data = StringIO()
    tree.export_graphviz(best_grid.estimators_[i-1], out_file=dot_data,feature_names=variable.columns)
    graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
    graph.write_png('./RF_result/tree' + str(i) + '.png')
print('Done')

まとめ

今回は、Random forestを使って疾患分類をしてみた。 我々の実験だと75.4%で先行研究の74.0%より1.4%高かった。 先行研究ではGray matter volumeのみであったが、我々はdiffusion MRIから得られるmapを入力。 とりあえず、SensitivityとSpecificityを調べたい。



頑張れ!喝!!の代わりにB!ブックマークを押していただけるとただただうれしいです(^^)! ↓