【Python初心者】入院日数を予測してみた
はじめに
AI は、入院時に患者の入院日数を予測するのに役立ちます。AIで入院期間を予測することでリスクの高い患者に的を絞った治療介入を行うことができる。また、医療従事者の負担軽減にもつながると考えられる。
病院に入院している人の入院期間を過去の入院データから予測してみました。MacBook Pro、Google Colaboratoryを使用しています。
１.必要なデータを手動で情報を取得
カルテから必要なデータをエクセルに入力する。それをCVSに変換する。
２.データの前処理
新しい特徴量の作成、余分な文字を削除、欠損値の処理などを行っています。
＃新しい特徴量の作成
def create_features(data: pd.DataFrame) -> pd.DataFrame:
data["挿管入院日一致"] = data["入院日"] == data["挿管日"]
# 入院日から転入日までの日数を作る
return datadata = create_features(raw_data) ＃学習データと訓練データに分ける
from sklearn.model_selection import train_test_split
X_raw = data.drop(["入院期間"], axis=1)
y = data["入院期間"]
_X_train, _X_test, y_train, y_test = train_test_split(X_raw, y, test_size=0.2)CategoryEndodersライブラリを使ってカテゴリカルな特徴量を作成しました。
参照 https://qiita.com/Hyperion13fleet/items/afa49a84bd5db65ffc31
import category_encoders as ce
from sklearn.preprocessing import RobustScaler
# データの前処理をする関数
def preprocess(data: pd.DataFrame, encoder = None, scaler = None, encode_type="one_hot") -> pd.DataFrame:
unused_cols = ["No.", "入院日", "生年月日", "Unnamed: 17", "Unnamed: 22", "主治医.1", "レベル", "4階転入日", "RST介入日", "RST終了日", "挿管日", "抜管日", "退院日"]
category_cols = ["性別", "主治医", "病名", "挿管", "肺炎", "挿管入院日一致", "吸引回数"]
# いらない列を削除
X = data.drop(unused_cols, axis=1)
# 欠損値の処理
X['人工呼吸器装着期間'] = X['人工呼吸器装着期間'].fillna(0)
# encoding
if encode_type == "one_hot":
if not encoder:
encoder = ce.one_hot.OneHotEncoder(cols=category_cols, use_cat_names=True)
encoder = encoder.fit(X)
X = encoder.transform(X)
elif encode_type == "label":
if not encoder:
encoder = ce.ordinal.OrdinalEncoder(cols=category_cols)
encoder = encoder.fit(X)
X = encoder.transform(X)
else:
print("ERROR: Unknown encode type")
print(X.info())
# scaling numeric features
if scaler == None:
scaler = RobustScaler()
scaler = scaler.fit(X)
X = scaler.transform(X)
X = pd.DataFrame(X, columns=scaler.get_feature_names_out())
return X, encoder, scalerX_train, encoder, scaler = preprocess(_X_train)
X_test, _, _ = preprocess(_X_test, encoder, scaler)
df_X, _, _= preprocess(X_raw)
df_y = y.to_numpy()
X_train_2, encoder, scaler = preprocess(_X_train)
X_test_2, _, _ = preprocess(_X_test, encoder, scaler, encode_type="label")
df_X_2, _, _= preprocess(X_raw, encode_type="label")
df_y_2 = y.to_numpy()３.取得した情報を元に機械学習モデルを作成
目的変数：入院期間
説明変数：性別,、主治医、病名、入院日、退院日等々
線形モデル
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_absolute_error,mean_squared_error
import numpy as np
model = LinearRegression()model.fit(X_train, y_train)#回帰
model = LinearRegression()
model.fit(df_X, df_y)
pred_y = model.predict(df_X)
print('決定係数（r2）:{}'.format(round(r2_score(df_y, pred_y),3)))
print('平均誤差（MAE）:{}'.format(round(mean_absolute_error(df_y, pred_y),3)))
print('RMSE:{}'.format(round(np.sqrt(mean_squared_error(df_y, pred_y)),3)))
重回帰分析（最小二乗法）
import statsmodels.api as sm
df_X = sm.add_constant(df_X)
model = sm.OLS(df_y, df_X)
result = model.fit()
print(result.summary()) OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.272
Model: OLS Adj. R-squared: 0.100
Method: Least Squares F-statistic: 1.577
Date: Tue, 10 May 2022 Prob (F-statistic): 0.0652
Time: 11:53:57 Log-Likelihood: -590.72
No. Observations: 121 AIC: 1229.
Df Residuals: 97 BIC: 1297.
Df Model: 23
Covariance Type: nonrobust
=================================================================================
coef std err t P>|t| [0.025 0.975]
---------------------------------------------------------------------------------
const 38.8647 9.161 4.242 0.000 20.683 57.047
年齢 -4.8895 4.788 -1.021 0.310 -14.392 4.613
性別_女性 -1.1533 3.834 -0.301 0.764 -8.762 6.456
性別_男性 1.1533 3.834 0.301 0.764 -6.456 8.762
主治医_I 11.0359 11.699 0.943 0.348 -12.183 34.255
主治医_H 11.0319 13.737 0.803 0.424 -16.233 38.297
主治医_F 11.4728 17.234 0.666 0.507 -22.731 45.677
主治医_E 6.7164 9.504 0.707 0.481 -12.146 25.579
主治医_C -8.1497 12.758 -0.639 0.524 -33.471 17.172
主治医_G 2.8025 8.522 0.329 0.743 -14.112 19.717
主治医_J -4.7695 7.147 -0.667 0.506 -18.954 9.415
主治医_A 5.3000 18.183 0.291 0.771 -30.788 41.388
主治医_B 3.4244 16.208 0.211 0.833 -28.743 35.592
病名_脳梗塞 16.4517 8.285 1.986 0.050 0.008 32.895
病名_くも膜下出血 -0.0909 13.053 -0.007 0.994 -25.998 25.816
病名_脳出血 8.8802 8.411 1.056 0.294 -7.813 25.573
病名_慢性硬膜下血腫 -1.5543 19.815 -0.078 0.938 -40.881 37.772
病名_頭部外傷 -5.6472 12.002 -0.471 0.639 -29.468 18.174
病名_脊椎疾患 17.8688 25.607 0.698 0.487 -32.955 68.692
病名_その他 2.9564 25.208 0.117 0.907 -47.075 52.987
挿管_有 -6.2833 16.293 -0.386 0.701 -38.621 26.054
挿管_無 -14.3231 17.647 -0.812 0.419 -49.348 20.702
挿管_nan 20.6064 32.493 0.634 0.527 -43.883 85.096
RST介入期間 8.6657 3.555 2.438 0.017 1.610 15.722
人工呼吸器装着期間 1.1040 2.620 0.421 0.674 -4.096 6.304
吸引回数_それ以外 -3.5072 3.915 -0.896 0.373 -11.277 4.262
吸引回数_多い 3.5072 3.915 0.896 0.373 -4.262 11.277
肺炎_無 -5.1125 4.190 -1.220 0.225 -13.429 3.204
肺炎_有 5.1125 4.190 1.220 0.225 -3.204 13.429
挿管入院日一致_True 0.3850 5.281 0.073 0.942 -10.097 10.867
挿管入院日一致_False -0.3850 5.281 -0.073 0.942 -10.867 10.097
==============================================================================
Omnibus: 89.754 Durbin-Watson: 2.053
Prob(Omnibus): 0.000 Jarque-Bera (JB): 696.470
Skew: 2.501 Prob(JB): 5.80e-152
Kurtosis: 13.635 Cond. No. 1.68e+16
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 2.78e-30. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.T値2以上。
p値は0.05以下なら、「有意性が高い」ですが、結果からはどの特徴量も有意性が高いものがありません。
統計的にRST介入期間だけ有意な結果が得られました。
４.作成したモデルを用いて入院日数を予測
import matplotlib.pyplot as plt
import japanize_matplotlib
#予測実測プロットの作成
plt.figure(figsize=(6,6))
plt.scatter(df_y, pred_y,color='blue',alpha=0.3)
y_max_ = max(df_y.max(), pred_y.max())
y_min_ = min(df_y.min(), pred_y.min())
y_max = y_max_ + (y_max_ - y_min_) * 0.1
y_min = y_min_ - (y_max_ - y_min_) * 0.1
plt.plot([y_min , y_max],[y_min, y_max], 'k-')
plt.ylim(y_min, y_max)
plt.xlim(y_min, y_max)
plt.xlabel('入院期間(observed)',fontsize=20)
plt.ylabel('入院期間(predicted)',fontsize=20)
plt.legend(loc='best',fontsize=15)
plt.title('yyplot(入院期間)',fontsize=20)
plt.savefig('yyplot.png')
plt.show()
print('回帰係数：',round(model.intercept_,3))
"""output
回帰係数： 11.072
"""
#回帰係数を格納したpandasDataFrameの表示
df_coef = pd.DataFrame({'coefficient':model.coef_.flatten()}, index=df_X.columns)
#グラフの作成
x_pos = np.arange(len(df_coef))
fig = plt.figure(figsize=(12,12))
ax1 = fig.add_subplot(1, 1, 1)
ax1.barh(x_pos, df_coef['coefficient'], color='b')
ax1.set_title('coefficient of variables',fontsize=18)
ax1.set_yticks(x_pos)
ax1.set_yticks(np.arange(-1,len(df_coef.index))+0.5, minor=True)
ax1.set_yticklabels(df_coef.index, fontsize=14)
ax1.set_xticks(np.arange(-20,20, 1))
ax1.set_xticklabels(np.arange(-20,20, 1),fontsize=12, rotation=90)
ax1.grid(which='minor',axis='y',color='black',linestyle='-', linewidth=1)
ax1.grid(which='major',axis='x',linestyle='--', linewidth=1)
plt.savefig('coef_sklearn.png',bbox_inches='tight')
plt.show()
別のモデルも試してみたランダムフォレスト（Label encoding）
# Random Forest (Label encoding)
from sklearn.ensemble import RandomForestRegressor
model_2 = RandomForestRegressor()
model_2.fit(X_train_2, y_train)
# trainデータに対しての予測とそれぞれのスコア
pred_y = model_2.predict(X_train_2)
print('決定係数（r2）:{}'.format(round(r2_score(y_train, pred_y),3)))
print('平均誤差（MAE）:{}'.format(round(mean_absolute_error(y_train, pred_y),3)))
print('RMSE:{}'.format(round(np.sqrt(mean_squared_error(y_train, pred_y)),3)))
# testデータに対しての予測とそれぞれのスコア
pred_y = model_2.predict(X_test_2)
print('決定係数（r2）:{}'.format(round(r2_score(y_test, pred_y),3)))
print('平均誤差（MAE）:{}'.format(round(mean_absolute_error(y_test, pred_y),3)))
print('RMSE:{}'.format(round(np.sqrt(mean_squared_error(y_test, pred_y)),3)))
結果を見ると明らかに過学習している結果になった。
今回のモデルで入院日数を予測することは難しそうなので、長期入院するリスクが何％あるかを予測するアプリを作成することにした。
新しい特徴量の生成
def create_features(data: pd.DataFrame) -> pd.DataFrame:
data["挿管入院日一致"] = data["入院日"] == data["挿管日"]
data["長期入院リスク"] = data["入院期間"] >= 60
data["人工呼吸器装着"] = data["人工呼吸器装着期間"].isnull()データの前処理
from sklearn.model_selection import train_test_split
X_raw = data.drop(["長期入院リスク"], axis=1)
y = data["長期入院リスク"]
_X_train, _X_test, y_train, y_test = train_test_split(X_raw, y, test_size=0.2)
import category_encoders as ce
from sklearn.preprocessing import RobustScaler
# データの前処理をする関数
def preprocess(data: pd.DataFrame, encoder = None, scaler = None, encode_type="one_hot") -> pd.DataFrame:
unused_cols = ["入院期間","人工呼吸器装着期間", "RST介入期間", "No.", "入院日", "生年月日", "Unnamed: 17", "Unnamed: 22", "主治医.1", "レベル", "4階転入日", "RST介入日", "RST終了日", "挿管日", "抜管日", "退院日"]
category_cols = ["性別", "主治医", "病名", "挿管", "肺炎", "挿管入院日一致", "吸引回数"]
# いらない列を削除
X = data.drop(unused_cols, axis=1)
# 欠損値補完
X["挿管"] = X["挿管"].fillna("無")
# encoding
if encode_type == "one_hot":
if not encoder:
encoder = ce.one_hot.OneHotEncoder(cols=category_cols, use_cat_names=True)
encoder = encoder.fit(X)
X = encoder.transform(X)
elif encode_type == "label":
if not encoder:
encoder = ce.ordinal.OrdinalEncoder(cols=category_cols)
encoder = encoder.fit(X)
X = encoder.transform(X)
else:
print("ERROR: Unknown encode type")
print(X.info())
# scaling numeric features
if scaler == None:
scaler = RobustScaler()
scaler = scaler.fit(X)
X = scaler.transform(X)
X = pd.DataFrame(X, columns=scaler.get_feature_names_out())
return X, encoder, scalerX_train, encoder, scaler = preprocess(_X_train)
X_test, _, _ = preprocess(_X_test, encoder, scaler)
df_X, _, _= preprocess(X_raw)
df_y = y.to_numpy()
X_train_2, encoder, scaler = preprocess(_X_train)
X_test_2, _, _ = preprocess(_X_test, encoder, scaler, encode_type="label")
df_X_2, _, _= preprocess(X_raw, encode_type="label")
df_y_2 = y.to_numpy()Logistic回帰モデル
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss, roc_auc_score
import numpy as np
model = LogisticRegression()
model.fit(X_train, y_train)
model.score(X_train, y_train)
0.7395833333333334model.score(X_test, y_test)
0.6#回帰
model = LogisticRegression()
model.fit(df_X, df_y)
pred_y = model.predict(df_X)
print('AUC:{}'.format(round(roc_auc_score(df_y, pred_y),3)))
print('Log loss:{}'.format(round(log_loss(df_y, pred_y),3)))
AUC:0.621
Log loss:8.849import pickle
# save the model to disk
filename = 'lr_model.pickle'
pickle.dump(model, open(filename, 'wb'))
# some time later...
# load the model from disk
loaded_model = pickle.load(open(filename, 'rb'))
result = loaded_model.score(df_X, df_y)
print(result)
0.743801652892562Random Forestモデル
# Random Forest (Label encoding)
from sklearn.ensemble import RandomForestClassifier
model_2 = RandomForestClassifier()
model_2.fit(X_train_2, y_train)# testデータに対しての予測とそれぞれのスコア
pred_y = model_2.predict(X_test_2)
print('AUC:{}'.format(round(roc_auc_score(y_test, pred_y),3)))
print('Log loss:{}'.format(round(log_loss(y_test, pred_y),3)))
AUC:0.583
Log loss:13.816今回、まず最初の目標は入院日数の予測を目標にモデル構築してきましたが、良い結果が得られず、長期入院する可能性に目標を変更しましたが、それでも結果は良くなかったです。やはり、データ数が少なかった事、入院に影響を与える因子の項目が少なかった事、新しい特徴量が少なかった事などがモデルの精度が上がらなかった原因だと考えられます。
５．アプリ作成
import os
import pandas as pd
from flask import Flask, request, render_template
import pickle
app = Flask(__name__)
# 学習済みモデルをロード
# save the model to disk
model_filename = 'lr_model.pickle'
encoder_filename = 'one_hot_encoder.pickle'
scaler_filename = 'scaler.pickle'
# load the model from disk
model = pickle.load(open("./model/" + model_filename, 'rb'))
encoder = pickle.load(open("./model/" + encoder_filename, 'rb'))
scaler = pickle.load(open("./model/" + scaler_filename, 'rb'))
@app.route('/', methods=['GET', 'POST'])
def form():
if request.method == 'POST':
data = pd.DataFrame(
{
"年齢": request.form['age'],
"性別": request.form['gender'],
"主治医": request.form['doctor'],
"病名": request.form['disease'],
"挿管": request.form['intubation'],
"吸引回数": request.form['suctions'],
"肺炎": request.form['pneumonia'],
"挿管入院日一致": request.form['intubation_hospitalization_day'],
"人工呼吸器装着": request.form['ventilator'],
},
index=[0],
)
X = encoder.transform(data)
X = scaler.transform(X)
print(X.T)
# 変換したデータをモデルに渡して予測する
predicted = model.predict_proba(X)[0][1] * 100
pred_answer = "長期入院になる確率は " + str(round(predicted, 2)) + "% です"
return render_template("index.html", answer=pred_answer)
return render_template("index.html", answer="")
if __name__ == "__main__":
port = int(os.environ.get('PORT', 8085))
app.run(host='0.0.0.0', port=port)Herokuにアプリをデプロイしました。
https://hosp-app-1.herokuapp.com
アプリは完成しましたが、精度はあまり良くありません。
このアプリの精度を上げるための課題として
データ数を増やす
一人の患者さんの説明変数を増やす
入院期間に影響がありそうな特徴量を増やす
どのようなデータが入院期間に影響するのか見つけるために多職種の人とコミュニケーションをとる
データ収集を続ける
6. おわりに
今回、プログラミングを学ぶのがほぼ初めてという状況からスタートしたのですが、実際にデータ使ってアプリを作成することは本当に良い経験ができました。まだまだ分からないことがたくさんありますが続けていきたいと思いました。
余談ですがこの間、たまたま大学病院に行く機会があり、そこでアンケートがありました。その内容はカルテをAIに学習させるためのデータ提供の協力のお願いでした。医療の分野でもAIが使われようとしているのだなと改めて実感しました。