患者の症状から入院期間を予測するモデルの作成

## はじめに
看護師として勤務していた際、時間的・人的リソースの不足により患者に提供できるケアが制限され、理想とする看護師像とのギャップに悩むことがありました。こうした現場での経験を通じて、AI・機械学習の技術を活用し、看護師の勤務環境に関する問題解決に取り組むことにしました。


## 解決したい社会問題
株式会社マイナビが2024/1/25～2024/4/17に実施した、2,608名の看護師および231件の事業者を対象にしたインターネット調査により、看護師の勤務環境とその問題が明らかになりました。

サービス残業・早出が「発生している」と回答した看護師は76.1％で前年比8.1ptの増加。時間外労働の主な業務は、「看護記録※」が57.7%で最多となり、次いで「患者対応・診療補助（56.1%）」、「患者情報収集・申し送り（34.0%）」となりました。
（※看護師が現場で実践した一連の看護ケアを記録したもの）
出典：https://www.mynavi.jp/news/2024/07/post_44358.html

本来、「看護記録」はケア実践後すぐに記録するべきですが、業務中は記録に充てる時間が確保できず、時間外で対応せざるを得ない現状があります。これは人員配置の不備や、過剰な業務量が主な原因です。さらに、シフト制であるにも関わらず「患者対応」までもが時間外に行われることは、業務の分担が適切でないことを示唆しています。


## 仮説
厚生労働省のデータによると看護職員数は年々増加していますが、看護師のサービス残業・早出や時間外労働は減少するどころか、むしろ増加しています。
![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3825323/72519894-7bf4-f99b-976e-b52d63f0b96f.png)
出典：https://www.mhlw.go.jp/content/10800000/001118192.pdf

この状況を打開するためには、単に人手を増やすだけでなく、業務の効率化と適切な人員配置が必要であることが分かります。
そこで、患者が入院する時点で、患者の症状に基づいて入院期間を予測できれば、今後のケア計画が立てやすくなり、それに伴って適切な人員配置が可能になるのではないかと考えました。これにより、看護師一人あたりの業務負担がより均等に分配され、結果的に時間外労働の削減につながるのではないかと考えました。


## 分析データ
「Hospital Length of Stay Dataset Microsoft」
Kaggleにて公開されている、入院患者の症状や入院期間に関する10万件のデータ。
出典：https://www.kaggle.com/datasets/aayushchou/hospital-length-of-stay-dataset-microsoft/data

データの内容を簡単に表でまとめると以下の通りです。
![スクリーンショット 2024-11-17 16.27.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3825323/2c654102-ef7d-6dad-0d2e-e58956d57385.png)


症状の項目に関しては、症状が無い場合は0、症状がある場合は1で表されています。
secondarydiagnosisnonicd9（ICD-9に該当しない二次診断）に関しては、WHOが勧告した９つに分かれた「国際疾病分類（ICD）」に該当しない症状を示しています。
血液データや、身体状態の項目に関しては、測定値が入力されています。
以上の情報を踏まえて、以下を取り決めてモデル作成をしていきます。

【目的】入院患者の状態から入院期間（lengthofstay）を予測する。
【目的変数】　lengthofstay
【説明変数】　eid、vdata、discharged、lengthofstayを除く全て。
【利用するアルゴリズム】線形回帰（LinearRegression）


## 実行環境
パソコン：MacBook Air
開発環境：Visual Studio Code 
言語：Python3.12.6
ライブラリ：Pandas、Numpy、Matplotlib、seaborn、Sklearn


## 分析の流れ
1.	データの確認
2.	データの前処理
3.	学習データとテストデータに分割
4.	予測モデルの学習
5.	予測モデルの評価


## 実行したコード

#### 1. データの確認 
まずは必要なライブラリのインポートと、データのCSVファイルを読み込みます。
```Python
# 必要なライブラリのインポート
!pip install numpy
!pip install pandas
!pip install seaborn
!pip install numpy
!pip install scikit-learn

import numpy as np 
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# データの読み込み
df = pd.read_csv('./data/LengthOfStay.csv')
puts 'code block.'
```

ydata_profilingを使用すればデータの分析結果をレポートとして出力してもらえますが、今回は勉強も兼ねて、自分でデータを確認していきます。
```Python
# データの最初の5行を表示
print(df.head())
# データの基本統計量を確認
df.describe()
# データの情報を確認
df.info()
```
lengthofstayの最小値は１日、最大値は１７日、平均値は４日だとわかり、目的変数の大まかな形をつかめました。タイプ型はint型、object型、float型の３種です。欠損値はありませんでした。
症状の項目がint型と判別されましたが、本来はbool型なので修正しておきます。
```Python
# 症状の列を抜き出す
symptoms = ['dialysisrenalendstage', 'asthma', 'irondef', 
           'pneum', 'substancedependence', 'psychologicaldisordermajor', 
           'depress', 'psychother', 'fibrosisandother',  'malnutrition',
           'hemo', 'secondarydiagnosisnonicd9']
# int型からbool型に変換
df[symptoms] = df[symptoms].astype(bool)
```

以上までの流れでデータの概要がつかめました。
更にこのデータの洞察を深めていきます。今回の目的変数であるlengthofstay と、説明変数との関係を確認しておこうと思います。（お急ぎの方は２まで飛んでください。）

インデックス、来院日、退院日は使用しなくても問題がなさそうなので、説明変数から外しておきます。
```Python
# 'eid', 'vdate', 'discharged'を除外する
df_filtered = df.drop(columns=['eid', 'vdate', 'discharged'])
print(df_filtered.columns)
```

まずは、基本情報のうち、カテゴリカルな列（gender、facid）とlengthofstayとの関係をみてみます。性別ごと、施設ごとでlengthofstayの平均値を棒グラフで可視化します。
```Python
# gender ごとの lengthofstay の平均を計算
gender_distinction = df_filtered.groupby('gender')['lengthofstay'].mean().reset_index()

# グラフの描画
fig, ax = plt.subplots(figsize=(6, 6))
sns.barplot(x='gender', y='lengthofstay', data=gender_distinction, ax=ax)
ax.set_title('Average lengthofstay by gender')
ax.set_xlabel('gender')
ax.set_ylabel('Average lengthofstay')

# 各バーの上に数値ラベルを追加して表示
for index, row in gender_distinction.iterrows():
    ax.text(index, row['lengthofstay'] + 0.05, round(row['lengthofstay'], 1), ha="center")
plt.show()

# facid ごとの lengthofstay の平均を計算
facid_distinction = df_filtered.groupby('facid')['lengthofstay'].mean().reset_index()

# グラフの描画
fig, ax = plt.subplots(figsize=(6, 6))
sns.barplot(x='facid', y='lengthofstay', data=facid_distinction, ax=ax)
ax.set_title('Average lengthofstay by facid')
ax.set_xlabel('facid')
ax.set_ylabel('Average lengthofstay')

# 各バーの上に数値ラベルを追加して表示
for index, row in facid_distinction.iterrows():
    ax.text(index, row['lengthofstay'] + 0.05, round(row['lengthofstay'], 1), ha="center")
plt.show()
```
![スクリーンショット 2024-11-17 15.31.46.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3825323/388aef8f-44dd-cd41-e83c-9bdcb427c66d.png)
genderとlengthofstayとの比較の結果、性差はほぼないことが分かりました。
facidとlengthofstayとの比較の結果を見ると、facid A,Bは入院期間が短く、facid C,D,Eの方が入院期間が長くなる傾向がありそうです。

facidごとのlengthofstayの平均値は分かりましたが、もしかすると外れ値や異常値が潜んでいるかもしれません。データの分布を見るためにfacidごとのヒストグラムを作成してみます。
```Python
# 2行3列のプロット領域を作成
fig, axs = plt.subplots(2, 3, sharex=True, sharey=True)
fig.suptitle('Histogram of lengthofstay by facid')

# facidごとのヒストグラムを作成
facilities = ['A', 'B', 'C', 'D', 'E']
for i, facid in enumerate(facilities):
    row, col = divmod(i, 3)  # facidごとにサブプロットの位置を決定
    axs[row, col].hist(df_filtered[df_filtered['facid'] == facid]['lengthofstay'], bins=20) # lengthofstayの最大値17を参考
    axs[row, col].set_title(f'facid {facid}')
    axs[row, col].set_xlabel('lengthofstay')
    axs[row, col].set_ylabel('Number of people')

# 不要なプロットを削除する
fig.delaxes(axs[1, 2])

# レイアウト調整して表示
plt.tight_layout()
plt.show()
```
![スクリーンショット 2024-11-17 15.35.01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3825323/b16d2d13-eb1e-09f5-f269-81e894756b51.png)

外れ値や異常値はありませんでした。facidCとDは母数が少ないことが分かりました。
やはりfacid AとBは入院期間が短い層が多いです。

つぎに、症状の列とlengthofstayとの関係をみてみます。
x軸をそれぞれの症状、y軸をlengthofstayの平均として棒グラフで可視化します。
```Python
# 各症状を持つ患者の平均入院日数を計算
symptom_distinction = []
for symptom in symptoms:
    avg_lengthofstay = df_filtered.loc[df_filtered[symptom] == 1, 'lengthofstay'].mean()
    symptom_distinction.append((symptom, avg_lengthofstay))

# データフレームに変換
symptom_distinction = pd.DataFrame(symptom_distinction, columns=['symptom', 'average_lengthofstay'])

# グラフの描画
plt.figure(figsize=(10, 6))
ax = sns.barplot(x='symptom', y='average_lengthofstay', data=symptom_distinction)
plt.title('Average lengthofstay by symptoms')
plt.xlabel('symptoms')
plt.ylabel('Average lengthofstay')
plt.xticks(rotation=90)

# 各バーの上に数値ラベルを追加
for index, row in symptom_distinction.iterrows():  
    ax.text(index, row['average_lengthofstay'] + 0.05, round(row['average_lengthofstay'], 1), ha="center")

# レイアウト調整して表示
plt.tight_layout()
plt.show()
```
![スクリーンショット 2024-11-17 15.37.16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3825323/ab658182-0fe0-2f12-6644-3958c83172d8.png)
dialysisrenalendstage（透析）やfibrosisandother（肺線維症）を有する患者は入院期間が長くなる傾向にありそうです。

最後に、数字が入力されている列とlengthofstayとの関係をみてみます。主に血液データと身体情報です。項目が多いため、相関行列のヒートマップを採用します。
その際、rcountの値の５＋はそのままだと扱いにくいため、５に置き換えます。
```Python
# 'rcount' 列で '5+' を 5 に変換
df_filtered['rcount'] = df_filtered['rcount'].replace('5+', 5)

# 数字列を抜き出す
numeric_columns = df_filtered[['lengthofstay', 'rcount', 'hematocrit', 'neutrophils', 'sodium', 'glucose', 'bloodureanitro', 'creatinine', 'bmi', 'pulse', 'respiration', 'secondarydiagnosisnonicd9']]

# 相関行列を作成
corr_matrix = numeric_columns.corr()

# ヒートマップ作成
plt.figure(figsize=(12, 8))
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', fmt=".2f")
plt.title('Correlation Heatmap of numeric_columns and lengthofstay')

# レイアウト調整して表示
plt.tight_layout()
plt.show()
```
![スクリーンショット 2024-11-17 15.39.47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3825323/91552718-2b6a-498b-373c-9346ba12a0e3.png)
結果、rcountが突出してlengthofstayとの相関が高いことがわかります。


#### 2. データの前処理
データへの理解が深まったところで、次はデータを分析しやすい形にしていきます。
カテゴリカルな列はそのままだと扱いにくいので、非数値データを数値化するためにOne-hotエンコーディングを行います。
```Python
# genderとfacidに対してOne-Hotエンコーディングを実装
encoded = pd.get_dummies(df_filtered[['gender', 'facid']], prefix=['gender', 'facid'])
# 元のデータフレームに結合
df_concat = pd.concat([df_filtered, encoded], axis=1)
# 'gender'と'facid'列を削除
df_concat.drop(['gender', 'facid'], axis=1, inplace=True)
# 結果を確認
df_concat.head()
```
これによりgenderはgender_F（女性）とgender_M（男性）に変換されました。
また、facidはfacid_ A〜Eの５つに変換されました。

つぎに、各要素のスケールを統一するために、標準化を行います。
本来、標準化は前処理にあたりますが、今回は学習データとテストデータに分割してから標準化を行いたいので、次の章で行います。


#### 3. 学習データとテストデータに分割 
今回の目的変数はlengthofstayなので、X（学習データ）はlengthofstay以外の項目、y（テストデータ）はlengthofstayとして、比率7:3で分割していきます。
```Python
# train_test_split関数をインポート
from sklearn.model_selection import train_test_split

# 学習データとテストデータに分割
X = df_concat.drop(columns=['lengthofstay']).to_numpy()
y = df_concat['lengthofstay'].to_numpy()

# 比率7:3で学習データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
print(X_train.shape, X_test.shape, y_train.shape, y_test.shape)
```
分割したデータX_trainとX_testのうち、bool型（症状）とOne-hotエンコーディングしたもの以外の項目を抽出して、標準化を行います。
X_testは、X_trainでfitさせたスケールを使用して変換していきます。
```Python
# StandardScalerクラスのインポート
from sklearn.preprocessing import StandardScaler

#　インデックスを確認
print(X_train.shape)
print(df_concat.drop(columns=['lengthofstay']).columns)

# 標準化したい説明変数を指定
columns_to_scale = [0,12,13,14,15,16,17,18,19,20,21]

# X_trainとX_testから、指定した説明変数を抽出
X_train_selected = X_train[:, columns_to_scale]
X_test_selected = X_test[:, columns_to_scale]

# StandardScalerを使って標準化
scaler = StandardScaler()
X_train_scaler = scaler.fit_transform(X_train_selected)
X_test_scaler = scaler.transform(X_test_selected)

# 元のデータをコピーし、標準化した部分だけを代入
X_train_scaled = np.copy(X_train)
X_train_scaled[:, columns_to_scale] = X_train_scaler
X_test_scaled = np.copy(X_test)
X_test_scaled[:, columns_to_scale] = X_test_scaler
```


#### 4. 予測モデルの学習 
1つの目的変数と、複数の説明変数を扱うので重回帰で分析していきます。
LinearRegressionクラスをインスタンス化してmodelという変数に代入し、標準化済みのX_train_scaledとX_test_scaledを利用して学習させます。
```Python
# 重回帰で分析
from sklearn.linear_model import LinearRegression
model = LinearRegression()
# 予測モデルの学習
model.fit(X_train_scaled, y_train)
```


 #### 5. 予測モデルの評価
まずは、決定係数を出力して予測モデルを評価していきます。
 ```Python
# 学習データを用いて決定係数を算出
model.score(X_train_scaled, y_train)
# テストデータを用いて決定係数を算出
model.score(X_test_scaled, y_test)
```
学習データは0.7614013358246663、テストデータは0.7634428403520678となりました。双方の決定係数に、大きな乖離は見られません。過学習も発生していないと判断できます。

次に、予測モデルの精度を総合的に評価するために、二乗平均平方根誤差（RMSE）と平均絶対誤差（MAPE）の2つの評価指標を用いて、誤差の絶対的な大きさと相対的な割合の両面からモデルの精度を確認します。
```Python
# テストデータに対する予測
y_pred = model.predict(X_test_scaled)

# 二乗平均平方根誤差（RMSE）の値を算出
from sklearn.metrics import mean_squared_error
mean_squared_error(y_test, y_pred, squared=False)

# 平均絶対誤差(MAPE) の値を算出
np.mean(np.abs((y_test - y_pred) / y_test)) * 100
```
二乗平均平方根誤差（RMSE）は1.1483359957786603となり、入院期間の誤差は約1.15日となりました。平均絶対誤差(MAPE)は32.09155514862062となり、妥当な予測精度と考えられます。

One-hotエンコーディングした項目を除外して、影響度も確認してみます。
```Python
# 計算式の確認
print("係数:", model.coef_)
print("切片:", model.intercept_)

# One-hotエンコーディングした項目以外のものを抽出
features = ['rcount', 'dialysisrenalendstage', 'asthma', 'irondef', 'pneum', 'substancedependence', 'psychologicaldisordermajor', 'depress', 'psychother', 'fibrosisandother', 'malnutrition', 'hemo', 'hematocrit', 'neutrophils', 'sodium', 'glucose', 'bloodureanitro', 'creatinine', 'bmi', 'pulse', 'respiration', 'secondarydiagnosisnonicd9']
coefficients = model.coef_[:22]

# 棒グラフのプロット
plt.figure(figsize=(10, 5))
ax = sns.barplot(x=features, y=coefficients)

# グラフの設定
plt.title('Correlation between lengthofstay and features')
plt.xlabel('features')
plt.ylabel('Correlation with Length of Stay')
plt.xticks(rotation=90)

# 各バーの上に数値ラベルを追加して表示
for index, value in enumerate(coefficients):
    ax.text(index, value + 0.03, round(value, 2), ha="center")
plt.show()
```
![スクリーンショット 2024-11-17 15.47.59.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3825323/0e0f5d01-94cc-a631-1822-2438ee2321f9.png)
ヒートマップでも相関が高かったrcountがやはり影響度が高そうです。
また、グラフ右側に位置している「血液データや身体状態などの項目」よりも、左側の「症状の項目」の方が影響度が高いことが分かりました。

最後に、実際の入院期間と、予測された入院期間の値をそれぞれ比較してみます。
```Python
# DataFrame に変換して表示
df = pd.DataFrame({
    'Actual': y_test,     # 実際の値
    'Predicted': y_pred   # 予測された値
})
df
```
![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3825323/8e2858b1-b93f-39a2-bc38-f31330e43a23.png)

以上の結果より、入院期間の予測モデルとして一定の有効性を確認することができました。


## 傾向と考察
facid（病院）ごとの入院期間の差は、病院規模の違いが主な要因だと考えられます。大きな病院は設備や環境が整っており、治療が長引く重症患者も受け入れることができるためです。
rcount（再入院回数）との相関が高い点に関しては、入退院を繰り返す病状であれば入院が長引くことが予想されるため、妥当な結果と言えます。
また、症状の項目は影響度が高く、dialysisrenalendstage（透析）、fibrosisandother（肺線維症）を有する患者は入院期間が長くなる傾向にあることがわかりました。
予測モデルの誤差は約1.15日であり、ケア計画やシフト作成にも十分に活用できると考えられます。

## 課題
データ対象者の国籍が不明であるため、この結果が日本人にも同様に適用できるかどうかは確かではありません。また、国によって治療方法が異なる場合があるため、その影響で入院期間に差が生じる可能性があります。


## まとめ
このような入院期間の予測モデルがあれば、ケア計画が立てやすくなり、業務の負荷を適切に分散できると思います。また、数値で可視化されることで忙しさへの心構えができ、精神的負担の軽減も期待できます。
今回の結果から、現場で実際に活用できる有意義なモデルが構築できたと考えています。現在のモデルで一定の成果が得られたため、さらなる精度向上を目指してLinearRegression 以外のアルゴリズムでの検証も実施していきたいです。


## 出典
株式会社マイナビ「看護師白書」
https://www.mynavi.jp/news/2024/07/post_44358.html
厚生労働省「看護職員就業者数の推移」
https://www.mhlw.go.jp/content/10800000/001118192.pdf
Kaggle「Hospital Length of Stay Dataset Microsoft」
https://www.kaggle.com/datasets/aayushchou/hospital-length-of-stay-dataset-microsoft/data