【Aidemy成果物】健診データによる肝疾患の判定-初学者による教師あり学習（分類）の実践-

YANA

はじめに

Aidemy「データ分析講座」の受講の成果物としてこの記事を作成します。

私は医療系の営業をしており、仕事のひとつとして患者のデータを収集して整理しては医師に渡すといった事をしておりました。
集めたデータでどんなことが言えるだろうかと考えることは好きでしたが、
データを統計学的に扱う知識がなかったので、データを分析して論文化するというレベルでは医師とディスカッションができず悶々としておりました。

そんなときにAidemyの存在を知り、医療系のデータを分析するスキルが身に付いたら今の仕事がさらに楽しくなりそうだな、それをメインの仕事にできたらさらに楽しそうだなと、ふわっとした期待を胸に講座の受講し今に至ります。

新しいスキルを身につけることは楽しいです。

概要

データとして肝疾患に関する健診データを用意しました。
私の仕事は心臓関係ですが、勉強の為に今回は専門外である肝臓領域のデータをSIGNATEから引っ張ってきました。

健診で得られたデータから肝疾患の有無を判断するモデルを作成します。

実行環境

Python 3.10.12
Google Colaboratory

分析の流れ

  1. データの読込

  2. データの確認

  3. データの整理

  4. モデルの作成と学習

コード

データの読込

使用するライブラリのインポートをします。

import numpy as np
import pandas as pd
import scipy
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set(font_scale=1)
from sklearn.model_selection import train_test_split
from sklearn.model_selection import KFold
from sklearn.model_selection import RandomizedSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from tqdm import tqdm
from sklearn.datasets import load_breast_cancer
from sklearn.metrics import accuracy_score
from sklearn.metrics import roc_curve
from sklearn.metrics import auc
!pip install optuna

import optuna

データの読込をします。
本記事内ではSIGNATEに提出するところまでコードを記述しないので評価用データ（test_df）,応募用サンプルファイル（sample_sub）は必要ないですが、自身のコンペ用の勉強も兼ねて一緒に扱っていきます。

train_df = pd.read_csv('/content/drive/MyDrive/Colab Notebooks/train.csv')
test_df = pd.read_csv('/content/drive/MyDrive/Colab Notebooks/test.csv')
sample_sub = pd.read_csv('/content/drive/MyDrive/Colab Notebooks/sample_submit.csv')

データの確認

読み込んだデータを確認していきます。

print(f'Train_df_shape : {train_df.shape}\n')
print(train_df.dtypes)
display(train_df.head())

学習用データは891例あり、Genderのみ数字ではないobject型ではない事がわかりました。

Train_df_shape : (891, 12)

id            int64
Age           int64
Gender       object
T_Bil       float64
D_Bil       float64
ALP         float64
ALT_GPT     float64
AST_GOT     float64
TP          float64
Alb         float64
AG_ratio    float64
disease       int64
dtype: object
インデックスが上から5番目までのデータ

数字の大小に意味はないのでIDをカテゴリカル変数を数値から文字列へ変更します。

train_df= train_df.astype({'id' : str})
test_df= test_df.astype({'id' : str})

数値データの統計量を表示して確認してみます。

display(train_df.describe())
display(test_df.describe())


train_dfの数値データ
test_dfの数値データ

train_dfのデータのAG_ratioにのみ欠損値がある事が確認できました。
普段扱っていない血液データの値なので数値自体にどのような意味があるのかは推察できませんが、T_Bil,D_Bil,ALP,ALT_GPTは平均値に対して最大値がかなり大きいので正規分布してなさそうと思います。
予想をしていく疾患の有無を表すdiseaseは平均値が0.5に近いのであまり偏りはないと判断します。

IDと性別のカテゴリカルデータの統計量も表示して確認します。

display(train_df.describe(exclude='number'))
display(test_df.describe(exclude='number'))


train_dfのカテゴリカルデータ
test_dfのカテゴリカルデータ


IDは全てのデータで重複していないことがわかりました。
データ内の男女差はありそうです。

訓練データとテストデータを確認できたので、これからは結合して各データを確認していきます。

all_df = pd.concat([train_df,test_df],axis=0).reset_index(drop=True)

all_df['Test_Flag'] = 0
all_df.loc[train_df.shape[0]: ,'Test_Flag'] = 1

欠損値があるAG_ratioを処理していきます。
欠損値を平均値で埋めることも考えたのですが、調べてみるとAG_ratioはAlb/TP-Albで計算されるそうなので（厳密に言うと違う）、欠損値はAlb/TP-Albで計算させて埋めます。
ちゃんと計算されてかつ他の値に変化が出ていないか確認したかったのでAG_ratio_newのカラムを作成

all_df['AG_ratio_new'] = all_df["AG_ratio"].mask(all_df["AG_ratio"].isna(), all_df['Alb']/(all_df['TP']-all_df['Alb']))

ちゃんと計算されたかの確認

na_row = all_df["AG_ratio"].isna()
all_df.loc[na_row, :]
AG_ratio_newの確認

AG_ratioに欠損がある行にAG_ratio_newとして計算された値が入りました。
心配なのでAG_ratioに欠損がないデータも一緒に確認します。

all_df.loc[[30,31,278,495,648,649],['TP','Alb','AG_ratio','AG_ratio_new']]


AG_ratioとAG_ratio_newが変わっていないか確認

AG_ratioはそのままAG_ratio_newに転記されたことが確認できました。

ここからグラフを使用してデータの確認と分析を行っていきます。

seabornを使用してデータを可視化します。
大抵の内科系の疾患は年齢が因子の一つになるので、データ内で年齢の偏りがないか確認します。

fig = sns.FacetGrid(all_df, col='disease', hue='disease', height=4)
fig.map(sns.histplot, 'Age', bins=30, kde=False)


疾患の有無に対する年齢分布

高齢者で特別多いわけではなく、疾患の有無で年齢の偏りはなさそうです。

各項目でヒストグラムを確認していきます。

fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(1, 5,figsize=(30,5))
sns.histplot(all_df['T_Bil'],kde=False,ax=ax1)
sns.histplot(all_df['D_Bil'],kde=False,ax=ax2)
sns.histplot(all_df['ALP'],kde=False,ax=ax3)
sns.histplot(all_df['ALT_GPT'],kde=False,ax=ax4)
sns.histplot(all_df['AST_GOT'],kde=False,ax=ax5)
血液データごとのヒストグラム①

データに偏りがあり外れ値が多そうなので箱ひげ図でも確認します。

fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(1, 5,figsize=(30,5))
sns.boxplot(y=all_df['T_Bil'],ax=ax1)
sns.boxplot(y=all_df['D_Bil'],ax=ax2)
sns.boxplot(y=all_df['ALP'],ax=ax3)
sns.boxplot(y=all_df['ALT_GPT'],ax=ax4)
sns.boxplot(y=all_df['AST_GOT'],ax=ax5)


血液データごとの箱ひげ図①

第一四分位と第三四分位の幅が狭く、第三四分位を超えるデータが多くあることが各項目で確認できました。

残りの項目も確認していきます。

fig, (ax1, ax2, ax3) = plt.subplots(1, 3,figsize=(30,5))
sns.histplot(all_df['TP'],kde=False,ax=ax1)
sns.histplot(all_df['Alb'],kde=False,ax=ax2)
sns.histplot(all_df['AG_ratio'],kde=False,ax=ax3)


血液データごとのヒストグラム②

こちらは正規分布してそうですが、箱ひげ図でも確認していきます。

血液データごとの箱ひげ図②

AG_ratioは第三四分位を超えるデータがある程度ありそうです。

次にヒートマップを使用して各項目の相関を確認していきます。

sns.heatmap(all_df[['disease','Age','T_Bil','D_Bil','ALP','ALT_GPT','AST_GOT','TP','Alb','AG_ratio_new']].corr(),
            vmax=1,vmin=-1,annot=True
            )


ヒートマップで各項目の相関を確認

ヒートマップからT_BilとD_Bil、ALT_GPTとAST_GOT、AlbがTPとAG_ratioに強い相関を示していることがわかりました。
相関しているカラムからどちらか１つ選んでモデルを作っていくことも考えましたが、せっかくデータがあるので工夫して全てのデータを活かしていきます。

データの整理

  • T_BilとD_Bil

細かいことは割愛しますが、T_Bil(総ビリルビン)はD_Bil(直接ビリルビン)とI_Bil(間接ビリルビン)の合計になり、疾患によりD_BilとI_Bilの増加量は異なるようです。そのためT_BilからD_Bilの値を引いてI_Bilを算出しD_BilとI_Bilを使用することにします。

  • ALT_GPTとAST_GOT

こちらも細かいことは割愛しますが、
AST、ALTはいずれも肝細胞の傷害の有無を推定するデータになり、ALTは主に肝臓に存在し、ASTは肝臓のみならず心筋や骨格筋、赤血球などにも広く存在するようです。そのため肝疾患があるとどちらも高値になります。ALTに加えてASTとALTの比を確認することで肝疾患の中で何の疾患かを判断する指標になるそうです。
以上からAST_GOTとALT_GPTのデータを用いてAST/ALTのデータを算出します。

  • TPとAlbとAG_ratio

TP(総たんぱく)はほとんどがAlb(アルブミン)とGlb(グロブリン)で構成されていて、健常者ではAlbが約67％、Glbが約33%で保たれており疾患によりバランスが崩れるそうです。このバランスがAG_ratioで表されAlb/Glbで算出されます。TPからAlbを引きGlbを算出し使用します。

all_df['I_Bil'] = all_df['T_Bil'] - all_df['D_Bil']
all_df['AST/ALT'] = all_df['AST_GOT'] / all_df['ALT_GPT']
all_df['Glb'] = all_df['TP'] - all_df['Alb']

新しく作成したカラムの相関を確認します。

sns.heatmap(all_df[['disease','Age','D_Bil','I_Bil','ALP','AST/ALT','AST_GOT','Glb','AG_ratio_new']].corr(),
            vmax=1,vmin=-1,annot=True
            )


新しい項目の相関を確認

各カラムの相関がかなり改善されました。
ここから学習に向けてコードを進めていきます。

モデルの作成と学習

今回は、ロジスティック回帰、RandomForest、LightGBMの3つのアルゴリズムを採用します。

  • ロジスティック回帰

外れ値が存在するとモデルの学習に悪影響を及ぼす可能性があるので、外れ値が多く確認されたカラムをカテゴリカル変数に変化させることにより、外れ値の影響を減らしていきます。

all_df['AgeBand'] = pd.qcut(all_df['Age'], 4)
all_df['D_BilBand'] = pd.qcut(all_df['D_Bil'], 4)
all_df['I_BilBand'] = pd.qcut(all_df['I_Bil'], 4)
all_df['ALPBand'] = pd.qcut(all_df['ALP'], 4)
all_df['AST_GOTBand'] = pd.qcut(all_df['AST_GOT'], 4)
all_df['AST/ALTBand'] = pd.qcut(all_df['AST/ALT'], 4)
all_df['AG_ratio_newBand'] = pd.cut(all_df['AG_ratio_new'], 4)

次にカテゴリカル変数に対しOne-Hot Encodingを施し、モデルが扱えるよう変換します。

all_df = pd.get_dummies(all_df, columns= ["Gender"])
all_df = pd.get_dummies(all_df, columns= ['AgeBand'])
all_df = pd.get_dummies(all_df, columns= ['D_BilBand'])
all_df = pd.get_dummies(all_df, columns= ['I_BilBand'])
all_df = pd.get_dummies(all_df, columns= ['ALPBand'])
all_df = pd.get_dummies(all_df, columns= ['AST_GOTBand'])
all_df = pd.get_dummies(all_df, columns= ['AST/ALTBand'])
all_df = pd.get_dummies(all_df, columns= ['AG_ratio_newBand'])

前処理を行うためにまとめていたall_dfのデータを訓練用とテスト用に再度分割します。

train = all_df[all_df['Test_Flag']==0]
test = all_df[all_df['Test_Flag']==1].reset_index(drop=True)

# 訓練データのdiseaseをtargetにする
target = train['disease']

# 学習に用いないカラムを削除
drop_col = ['id','Age','Test_Flag','disease','D_Bil','T_Bil','I_Bil','ALP','ALT_GPT','AST_GOT','TP','Alb','AG_ratio','AG_ratio_new']

train = train.drop(drop_col, axis=1)
test = test.drop(drop_col, axis=1)

学習に用いるデータを一通り準備できたのでロジスティック回帰を用いて予測モデルを構築していきます。
モデルの評価は今回KFold法を使用します。
また評価はaccuracyとAUCで評価します。

cv = KFold(n_splits=10, random_state=0, shuffle=True)

train_acc_list = []
val_acc_list = []
val_auc_list = []

#fold毎に学習データのインデックスと評価データのインデックスを得る

for i ,(trn_index, val_index) in enumerate(cv.split(train, target)):

  print(f'Fold ; {i}')

  X_train, X_val = train.loc[trn_index], train.loc[val_index]
  y_train, y_val = target[trn_index], target[val_index]

  model = LogisticRegression()
  model.fit(X_train, y_train)

  y_pred = model.predict(X_train)
  train_acc = accuracy_score(y_train, y_pred)

  print(train_acc)
  train_acc_list.append(train_acc)

  y_pred_val = model.predict(X_val)
  val_acc = accuracy_score(y_val, y_pred_val)

  print(val_acc)
  val_acc_list.append(val_acc)

  y_val_proba = model.predict_proba(X_val)
  pos_prob_val = model.predict_proba(X_val)[:, 1]
  fpr_val, tpr_val, thresholds = roc_curve(y_val, pos_prob_val)
  print(auc(fpr_val, tpr_val))
  val_auc_list.append(auc(fpr_val, tpr_val))


print('-'*10 + 'Result' + '-'*10)
print(f'Train_acc : {train_acc_list} , Ave : {np.mean(train_acc_list)}')
print(f'Valid_acc : {val_acc_list} , Ave : {np.mean(val_acc_list)}')
print(f'Valid_auc : {val_auc_list} , Ave : {np.mean(val_auc_list)}')
----------Result----------
Train_acc : [0.7940074906367042, 0.7992518703241895, 0.8104738154613467, 0.8017456359102244, 0.7955112219451371, 0.8029925187032418, 0.8054862842892768, 0.7942643391521197, 0.7892768079800498, 0.8104738154613467] , Ave : 0.8003483799863638
Valid_acc : [0.8, 0.8202247191011236, 0.6966292134831461, 0.7865168539325843, 0.8202247191011236, 0.8202247191011236, 0.7528089887640449, 0.8089887640449438, 0.8426966292134831, 0.7640449438202247] , Ave : 0.7912359550561798
Valid_auc : [0.891358024691358, 0.8863636363636364, 0.7847011144883486, 0.8465045592705166, 0.9146341463414633, 0.9082051282051282, 0.8292053663570691, 0.9035714285714286, 0.9155712841253791, 0.8348813209494325] , Ave : 0.871499600936376
y_val_proba = model.predict_proba(X_val)

# 陽性の確率だけが必要なので[:, 1]をして陰性の確率を落とす
pos_prob_val = model.predict_proba(X_val)[:, 1]
fpr_val, tpr_val, thresholds = roc_curve(y_val, pos_prob_val)
# 描画
plt.plot(fpr_val,tpr_val)
plt.xlabel('1-specificity (FPR)')
plt.ylabel('sensitivity (TPR)')
plt.title('ROC Curve')
plt.show()


ロジスティック回帰でのROC Curve

ロジスティック回帰では
ACC : 0.791, AUC : 0.871
でした。
次はRandomForestで予測モデルを作成していきます。

  • RandomForest

途中までロジスティック回帰と同じ処理をしていきます。

train_df = pd.read_csv('/content/drive/MyDrive/Colab Notebooks/train.csv')
test_df = pd.read_csv('/content/drive/MyDrive/Colab Notebooks/test.csv')
sample_sub = pd.read_csv('/content/drive/MyDrive/Colab Notebooks/sample_submit.csv')


#訓練データとテストデータを結合

all_df = pd.concat([train_df,test_df],axis=0).reset_index(drop=True)

all_df['Test_Flag'] = 0
all_df.loc[train_df.shape[0]: ,'Test_Flag'] = 1

#AG_ratioの欠損値に対して計算式で代入する

all_df['AG_ratio_new'] = all_df["AG_ratio"].mask(all_df["AG_ratio"].isna(), all_df['Alb']/(all_df['TP']-all_df['Alb']))

#相関している説明変数を工夫

all_df['I_Bil'] = all_df['T_Bil'] - all_df['D_Bil']
all_df['AST/ALT'] = all_df['AST_GOT'] / all_df['ALT_GPT']
all_df['Glb'] = all_df['TP'] - all_df['Alb']

RandomForestとLightGBMでは外れ値を多く含むカラムをカテゴリカル変数で扱うのではなく、対数変換して扱います。
決定木をベースとしたモデルは特徴量の大小関係のみに着目しており、値自体には意味がないので対数変換による正規化は必要ないかもしれませんが、自身の勉強用と対数変換した方がRandomForestは値が少し良かったので今回は対数変換して進めていきます。
AgeはlightGBMで扱いやすいようにqcutではなく、cut関数を使用しました。
GenderとAgeBandは同じようにOne-Hot_Encordingで変換します。

all_df['AgeBand'] = pd.cut(all_df['Age'], [0,20,40,60,80,100])

all_df['D_Bil_log'] = np.log10(all_df['D_Bil'] + 1)
all_df['I_Bil_log'] = np.log10(all_df['I_Bil'] + 1)
all_df['ALP_log'] = np.log10(all_df['ALP'] + 1)
all_df['AST_GOT_log'] = np.log10(all_df['AST_GOT'] + 1)
all_df['AST/ALT_log'] = np.log10(all_df['AST/ALT'] + 1)
all_df['AG_ratio_new_log'] = np.log10(all_df['AG_ratio_new'] + 1)

#Gender,AgeBandをOne-Hot_Encordingで変換

all_df = pd.get_dummies(all_df, columns= ["Gender"])
all_df = pd.get_dummies(all_df, columns= ["AgeBand"])
#対数変換したカラムをヒストグラムで確認
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4,figsize=(30,5))
sns.histplot(all_df['D_Bil_log'],kde=False,ax=ax1)
sns.histplot(all_df['ALP_log'],kde=False,ax=ax2)
sns.histplot(all_df['AST_GOT_log'],kde=False,ax=ax3)
sns.histplot(all_df['AG_ratio_new_log'],kde=False,ax=ax4)


対数変換後のヒストグラム


対数変換前のヒストグラム

all_dfのデータを訓練用とテスト用に再度分割します。

train = all_df[all_df['Test_Flag']==0]
test = all_df[all_df['Test_Flag']==1].reset_index(drop=True)

target = train['disease']

drop_col = ['id','Age','Test_Flag','disease','D_Bil','T_Bil','I_Bil','ALP','ALT_GPT','AST_GOT','TP','Alb','AG_ratio','AG_ratio_new']

train = train.drop(drop_col, axis=1)
test = test.drop(drop_col, axis=1)

RandomForestをランダムサーチを使用してパラメーターチューニングをします。グリッドサーチは時間がかかり過ぎたため途中で諦めました。
また試行錯誤しながら何回もコードを実行するにはKFoldだと時間がかかるためパラメーターチューニングはホールドアウト法で行いました。

#訓練データと検証データをホールドアウト法で作成
X_train, X_val, y_train, y_val = train_test_split(train,target,test_size = 0.2, random_state = 0)

# サーチするパラメータを設定
RF_params = {
    RandomForestClassifier(): {"n_estimators": scipy.stats.randint(10, 300),
                               "max_depth": scipy.stats.randint(1, 100),}
}

# スコア比較用に変数を設定
max_score = 0
best_param = None

# ランダムサーチでパラメーターサーチ
for model, param in RF_params.items():
    clf = RandomizedSearchCV(model, param)
    clf.fit(X_train, y_train)
    y_pred_val = clf.predict(X_val)
    score = accuracy_score(y_val, y_pred_val)
    
    if max_score < score:
        max_score = score
        best_param = clf.best_params_

print("Best_params:{}".format(best_param))
print("Best_score:",max_score)
Best_params:{'max_depth': 69, 'n_estimators': 236}
Best_score: 0.8547486033519553

得られたパラメーターを用いて、KFold法でも検証していきます。

#訓練データと検証データをKFold法で作成
cv = KFold(n_splits=10, random_state=0, shuffle=True)

train_acc_list = []
val_acc_list = []
val_auc_list = []

#fold毎に学習データのインデックスと評価データのインデックスを得る

for i ,(trn_index, val_index) in enumerate(cv.split(train, target)):

  print(f'Fold ; {i}')

  X_train, X_val = train.loc[trn_index], train.loc[val_index]
  y_train, y_val = target[trn_index], target[val_index]

  model_RF = RandomForestClassifier(random_state=0,max_depth=69,n_estimators=236)
  model_RF.fit(X_train, y_train)

  y_pred = model_RF.predict(X_train)
  train_acc = accuracy_score(y_train, y_pred)

  print(train_acc)
  train_acc_list.append(train_acc)

  y_pred_val = model_RF.predict(X_val)
  val_acc = accuracy_score(y_val, y_pred_val)

  print(val_acc)
  val_acc_list.append(val_acc)

  y_val_proba = model_RF.predict_proba(X_val)
  pos_prob_val = model_RF.predict_proba(X_val)[:, 1]
  fpr_val, tpr_val, thresholds = roc_curve(y_val, pos_prob_val)
  print(auc(fpr_val, tpr_val))
  val_auc_list.append(auc(fpr_val, tpr_val))


print('-'*10 + 'Result' + '-'*10)
print(f'Train_acc : {train_acc_list} , Ave : {np.mean(train_acc_list)}')
print(f'Valid_acc : {val_acc_list} , Ave : {np.mean(val_acc_list)}')
print(f'Valid_auc : {val_auc_list} , Ave : {np.mean(val_auc_list)}')
----------Result----------
Train_acc : [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] , Ave : 1.0
Valid_acc : [0.8555555555555555, 0.8651685393258427, 0.797752808988764, 0.8426966292134831, 0.8651685393258427, 0.9325842696629213, 0.898876404494382, 0.8764044943820225, 0.9325842696629213, 0.8876404494382022] , Ave : 0.8754431960049937
Valid_auc : [0.9375308641975308, 0.956060606060606, 0.940982776089159, 0.9354103343465046, 0.9596036585365854, 0.9838461538461537, 0.9674922600619195, 0.9704081632653062, 0.9656218402426695, 0.9385964912280702] , Ave : 0.9555553147874505

訓練データではぴったり予想が可能になり、
検証データでは
ACC : 0.875 AUC : 0.955
でした。
ROCカーブもグラフ化します。

pos_prob_val = model_RF.predict_proba(X_val)[:, 1]
fpr_val, tpr_val, thresholds = roc_curve(y_val, pos_prob_val)

plt.plot(fpr_val,tpr_val)
plt.xlabel('1-specificity (FPR)')
plt.ylabel('sensitivity (TPR)')
plt.title('ROC Curve')
plt.show()


RandomForest ROC Curve

次はlightGBMで予測モデルを作成していきます。
lightGBMはハイパーパラメータ探索ライブラリの「optuna」を用いて、パラメータを探索していきます。

def objective(trial):

    # RandomForestと同様にKFoldでモデルを評価
    cv = KFold(n_splits=10, random_state=0, shuffle=True)
    lgb_val_acc_list = []

    # ハイパーパラメータを定義
    lgb_params = {
        "objective":"binary",
        "metric": "binary_error",
        "force_row_wise" : True,
        "seed" : 0,
        'min_data_in_leaf': 5,

        # optunaでパラメータ探索する区間を指定
        "max_depth": trial.suggest_int('max_depth', 3, 20),
        "num_leaves": trial.suggest_int('num_leaves', 21, 41),
        "learning_rate": trial.suggest_uniform('learning_rate', 0.1, 0.2),
        }

    # LightGBMがカラム名に[]や{}を含んでいるとエラーが出るので、変更する
    rename_dict = {
        'AgeBand_(0, 20]' : 'AgeBand_1',
        'AgeBand_(20, 40]' : 'AgeBand_2',
        'AgeBand_(40, 60]' : 'AgeBand_3',
        'AgeBand_(60, 80]' : 'AgeBand_4',
        'AgeBand_(80, 100]' : 'AgeBand_5',
        }

    # KFold で学習させる
    for i ,(trn_index, val_index) in enumerate(cv.split(train, target)):

        print(f'Fold : {i}')
        X_train ,X_val = train.loc[trn_index].rename( columns =rename_dict), train.loc[val_index].rename( columns =rename_dict)
        y_train ,y_val = target[trn_index], target[val_index]

        lgb_train = lgb.Dataset(X_train, y_train)
        lgb_valid = lgb.Dataset(X_val, y_val)

        model_lgb = lgb.train(
            params = lgb_params,
            train_set = lgb_train,
            valid_sets = [lgb_train, lgb_valid],
            verbose_eval = -1 ,
            early_stopping_rounds=10
           )

        y_pred = model_lgb.predict(X_train)

        y_pred_val = model_lgb.predict(X_val)
        val_acc = accuracy_score(y_val, np.where(y_pred_val>=0.5, 1, 0))

        lgb_val_acc_list.append(val_acc)

    # 最適化の指標にする値を返り値に設定
    return np.mean(val_acc_list)

# 最適化タスクを定義
study = optuna.create_study(direction="maximize")

# 最適化を実行
study.optimize(objective, n_trials=100)

print("best_value", study.best_value)
print("best_params", study.best_params)
best_value 0.8810611735330836
best_params {'max_depth': 12, 'num_leaves': 21, 'learning_rate': 0.17172833900653073}

lightGBMが一番良いスコアになりました。

同様のコードでAUCも求めました。

def objective(trial):

    cv = KFold(n_splits=10, random_state=0, shuffle=True)
    lgb_val_auc_list = []

    lgb_params = {
        "objective":"binary",
        "metric": "binary_error",
        "force_row_wise" : True,
        "seed" : 0,
        'min_data_in_leaf': 5,

        # optunaでパラメータ探索する区間を指定
        "max_depth": trial.suggest_int('max_depth', 3, 20),
        "num_leaves": trial.suggest_int('num_leaves', 21, 41),
        "learning_rate": trial.suggest_uniform('learning_rate', 0.1, 0.2),
        }

        'AgeBand_(0, 20]' : 'AgeBand_1',
        'AgeBand_(20, 40]' : 'AgeBand_2',
        'AgeBand_(40, 60]' : 'AgeBand_3',
        'AgeBand_(60, 80]' : 'AgeBand_4',
        'AgeBand_(80, 100]' : 'AgeBand_5',
        }

    for i ,(trn_index, val_index) in enumerate(cv.split(train, target)):

        print(f'Fold : {i}')
        X_train ,X_val = train.loc[trn_index].rename( columns =rename_dict), train.loc[val_index].rename( columns =rename_dict)
        y_train ,y_val = target[trn_index], target[val_index]

        lgb_train = lgb.Dataset(X_train, y_train)
        lgb_valid = lgb.Dataset(X_val, y_val)

        model_lgb = lgb.train(
            params = lgb_params,
            train_set = lgb_train,
            valid_sets = [lgb_train, lgb_valid],
            verbose_eval = -1 ,
            early_stopping_rounds=10
           )

        y_pred = model_lgb.predict(X_train)

        fpr, tpr, thresholds = metrics.roc_curve(y_val, y_pred_val, pos_label=2)
        val_auc = metrics.auc(fpr, tpr)

        lgb_val_auc_list.append(val_auc)

    return np.mean(val_auc_list)

study = optuna.create_study(direction="maximize")

study.optimize(objective, n_trials=100)

print("best_value", study.best_value)
print("best_params", study.best_params)
best_value 0.9555553147874505
best_params {'max_depth': 4, 'num_leaves': 26, 'learning_rate': 0.12946386534426424}

AUCはRandomForestと同じスコアになりました。

考察・まとめ

ロジスティック回帰、RandomForest、lightGBMの３つのモデルを使用して健診データから肝疾患の予想を行いました。
だいぶ試行錯誤しましたが前処理を行うことで、予測精度が高まることが身をもって経験することができました。
一方コードにおいては、各モデルの作成、ハイパーパラメータのサーチ、各スコアの算出までの流れの部分がコードを細切れに区切ってなんとか形にしたという感じなので、不格好に感じます。
憧れのきれいなコードをパパっと書く姿になるためにはまだまだ時間がかかりそうです。

実践で試行錯誤しながらデータ分析に取り組むことで学習したことの理解が大幅に進むことが分かったので、これからもSIGNATEやKaggleの問題に果敢に挑戦していこうと思います。