はじめに

日本の実献血者数は年々減少傾向にあり、献血可能人口に占める割合は約6%程度¹にとどまっています。このような状況が続くと、輸血用の血液在庫確保に支障が生じ、需要に応じた血液在庫の確保が困難になると推測されます。献血事業の将来を見据えた献血シミュレーション²やマルコフモデル³を用いた予測が行われていることから予測のニーズは高いと考えました。しかし、月次レベルでの予測を行う研究を見つけることができなかったので、本記事では、季節変動自己回帰移動平均（SARIMA）モデルを用いて月次レベルでの予測モデルの開発を試みます。

予測結果は以下のサイトで公開しています。

献血者のデータ

使用するすべてのデータは日本赤十字社の公式ホームページに掲載されている献血者数速報⁴より2017年1月から2024年4月（執筆時点）までのデータを取得しています。

月ごとの全国の総献血者数をデータセットとして用意しました。

ここから、献血者のデータを可視化・分析をしていきます。

変動成分

# トレンド、季節性、残差を求める（加法モデル）
res_additive = sm.tsa.seasonal_decompose(df.values, period=12, model='additive')

# グラフ化
fig = res_additive.plot()
plt.show()

原形列に対して変動成分の分解を行いました。一番上のグラフから生データ、データの長期的な傾向を表すトレンド、データの周期的な変動を表す季節変動、トレンドと季節変動でとらえることができない成分である残差をプロットしています。

この結果から季節変動に注目して、どのくらいの周期で増減を繰り返しているのかを確認します。

s = res_additive.seasonal[60]

additive_seasona = res_additive.seasonal[60:84]

# トレンド、季節性、残差を求める（加法モデル）
plt.plot(df.index[60:84], additive_seasona)
plt.axhline(s, ls = "-.", color = "magenta")
plt.show()

試しに2022年1月から2024年1月までのデータで確認してみると、12か月周期で繰り返されていそうです。

定常性

定常性とは時間によらず統計的な特性（平均、分散、自己相関など）が一定である性質のことです。ここでは原形列に対して1期前のデータとの差分を計算し階差系列を得る差分変換を行います。

# 差分変換
df_diff = grouped_date_total.diff(1).dropna()
plt.plot(df_diff)
plt.show()

ここで、定常性を確認する手法として差分系列に対してADF検定を行います。

ADF検定の流れ：

帰無仮説を対立仮説を立てる
- 帰無仮説$H_0$：時系列は定常である（単位根を持つ）
- 対立仮説$H_1$：時系列は定常ではない（単位根を持たない）
検定量の比較

一般的にp値の閾値(有意水準)として0.05が用いられます。p値が0.05より小さいとき5%有意といい、帰無仮説を棄却し対立仮説を採択します。また、ADF統計量が5%の棄却限界値より小さいとき5%有意といい、帰無仮説を棄却し対立仮説を採択します。

すなわち、p値とADF統計量の両方が5%の有意水準を下回る場合、時系列は定常である可能性が高いことを示唆します。

# ADF検定
dftest = adfuller(df_diff)
print("--- diff ---")
print('ADF Statistic: %f' % dftest[0])
print('p-value: %f' % dftest[1])
print('Critical values :')
for k, v in dftest[4].items():
    print('\t', k, v)

--- diff ---
ADF Statistic: -4.775412
p-value: 0.000061
Critical values :
	 1% -3.5194805351545413
	 5% -2.9003945086747343
	 10% -2.5874984279778395

ADF統計量は-4.775412であり5%の棄却限界値を下回っているので5%有意
p値は0.000061であり5%の棄却限界値を下回っているので5%有意

よって、帰無仮説（時系列は非定常である）は棄却され、対立仮説（時系列は定常である）が採択されました。すなわち、定常性が認められたということです。👏

相関分析

相関分析とは、2つ以上の変数間の関係性を定量的に評価するための統計手法で、変数間の関連性の強さと方向性を理解することができます。
自己相関（ACF）：
さまざまなラグ（遅れ）での自己相関を示す関数。ACFプロットを使用することで、時系列データ内の観測値が過去の自身とどの程度関連しているかを視覚的に評価することができます。
偏自己相関（PACF）：
中間のラグの影響を除いた特定のラグでの偏自己相関を示す関数。偏自己相関関数を用いることで、直接的なラグの影響を評価することができます。

# サブプロットの作成
fig, axes = plt.subplots(1, 2)

# 差分変換したデータのACFとPACFをプロット
plot_acf(df_diff, lags=36, ax=axes[0])
plot_pacf(df_diff, lags=36, ax=axes[1])

# グラフを表示
axes[0].set_title('ACF')
axes[1].set_title('PACF')
plt.tight_layout()
plt.show()

これらのコレログラムより、青い部分を帰無仮説（相関がない）とした 95% 信頼区間を表しています。つまり、青い部分の外にある値は 5% の有意水準で相関があることになります。よって、自己相関(ACF)ではラグ5,7,12で相関があることが確認できます。

モデルの構築と予測

SARIMAモデル

SARIMAモデルは一般的に以下の式(1)で表すことができます。

$$
SARIMA(p,d,q)(P,D,Q)_s\qquad(1)
$$

(P)p-(季節性)自己回帰項の次数：

自己回帰(AR)のパラメータ。現在の値が過去の値のどれだけに依存するかを示す。

(D)d-(季節性)階差の次数：

和分(I)のパラメータ。時系列データを定常にするために差分を取る回数を示す。

(Q)q-(季節性)移動平均項の次数：

移動平均(MA)のパラメータ。過去の誤差が将来の値にどれだけ影響を与えるかを示す。

s-季節変動の周期：

季節性の周期を示す。

AIC

式(1)のそれぞれの次数に対し0か1を考えるとしても、 2^7=128 通りのモデルを考えなければなりません。そのため、周期sは自己相関関数をもとに決めたり、赤池情報量基準(AIC)を用いて、モデルの最尤推定し、AICが最も小さくなるようなモデルを選択するといった方針を立てます。

AICは一般的に以下の式(2)で表すことができます。

$$
AIC=-2\log{L}+2K\qquad(2)
$$

L-モデルの最大尤度：

モデルがデータにどれだけ適合しているかを示す指標⁵

K-パラメータ数：

説明変数の数

AICが小さいほど良いモデルといえますが、値そのものに意味はありません。複数のモデルを比較したとき、AICの小さいほうがより良いモデルという相対的な評価をします。⁶

パラメータチューニング

自己相関関数より周期sは12、定常性検定より階差の次数d(D)は1としておきます。

ここで、auto_arima関数を用いて最適パラメータの探索を行います。

# SARIMAモデルの自動パラメータ探索
smodel = pm.auto_arima(df_train, start_p=0, start_q=0, max_p=3, max_q=3, 
                       seasonal=True, m=12, start_P=0, start_Q=0, 
                       max_P=3, max_Q=3, d=1, D=1, trace=True,
                       error_action='ignore', suppress_warnings=True, 
                       stepwise=True)

Performing stepwise search to minimize aic
 ARIMA(0,1,0)(0,1,0)[12]             : AIC=1406.214, Time=0.01 sec
 ARIMA(1,1,0)(1,1,0)[12]             : AIC=1408.423, Time=0.07 sec
 ARIMA(0,1,1)(0,1,1)[12]             : AIC=1412.075, Time=0.09 sec
 ARIMA(0,1,0)(1,1,0)[12]             : AIC=1406.478, Time=0.03 sec
 ARIMA(0,1,0)(0,1,1)[12]             : AIC=1406.515, Time=0.04 sec
 ARIMA(0,1,0)(1,1,1)[12]             : AIC=inf, Time=0.46 sec
 ARIMA(1,1,0)(0,1,0)[12]             : AIC=1407.906, Time=0.01 sec
 ARIMA(0,1,1)(0,1,0)[12]             : AIC=1412.027, Time=0.04 sec
 ARIMA(1,1,1)(0,1,0)[12]             : AIC=1414.185, Time=0.07 sec
 ARIMA(0,1,0)(0,1,0)[12] intercept   : AIC=1407.514, Time=0.01 sec

Best model:  ARIMA(0,1,0)(0,1,0)[12]          
Total fit time: 0.838 seconds

よって、$SARIMA(p,d,q)(P,D,Q)_s=SARIMA(0,1,0)(0,1,0)12$と求まりました。

結果

テストデータに対して構築したSARIMAモデルで予測を行います。

# 学習
sarima_model = SARIMAX(train, order=(0, 1, 0), seasonal_order=(0, 1, 0, 12))
sarima_fit = sarima_model.fit()

# 予測
#学習データの期間の予測値
train_pred = sarima_fit.predict()
#テストデータの期間の予測値
test_pred = sarima_fit.forecast(len(df_test))

実観測値の学習データを青色、テストデータを灰色でプロットし、SARIMAモデルの予測値を水色でプロットさせます。

# グラフ化
fig, ax = plt.subplots()
ax.plot(df_train[24:].index, df_train[24:].values, label="actual(train dataset)")
ax.plot(df_test.index, df_test.values, label="actual(test dataset)", color="gray")
ax.plot(df_train[24:].index, train_pred[24:].values, color="c")
ax.plot(df_test.index, test_pred.values, label="SARIMA", color="c") 
ax.legend()
plt.show()

増加・減少傾向はある程度つかめているものの、大きく予測が外れてしまっている月があります。

精度評価

ここまでは、グラフを見て定性的に精度の評価をしましたが、定量的に精度評価を実施します。

時系列予測や回帰問題においてモデルの予測性能を評価するための評価関数としてしばしば用いられるRMSE、MAE、およびMAPEを計算します。

RMSE (Root Mean Squared Error, 平均二乗平方根誤差)

一般的に以下の式(3)で定義されます。⁷

$$
\text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2}\qquad(3)
$$

誤差の二乗平均の平方根を取ったもので、予測値と実測値の間の差を示します。誤差の大きさに比例する形で影響を与えるため、大きな誤差が強調されます。
MAE (Mean Absolute Error, 平均絶対誤差)

一般的に以下の式(4)で定義されます。⁷

$$
\text{MAE} = \frac{1}{n} \sum_{i=1}^{n} \left| y_i - \hat{y}_i \right|\qquad(4)
$$

誤差の絶対値の平均を取ったもので、予測値と実測値の間の平均的な差を示します。外れ値の影響を受けにくい指標です。
MAPE (Mean Absolute Percentage Error, 平均絶対パーセント誤差)

一般的に以下の式(5)で定義されます。⁷

$$
\text{MAPE} = \frac{100}{n} \sum_{i=1}^{n} \left| \frac{\hat{y}_i - y_i}{y_i} \right|\qquad(5)
$$

実測値に対する誤差の割合を示します。誤差が実測値に対してどれくらいの割合であるかを示すため、スケールに依存しない指標です。

このとき、

$n$：データ数

${y}_i$：i番目の実測値

$\hat{y}_i$：i番目の予測値

を表します。

ここではそれぞれscikit-learnのmean_squared_error関数、mean_absolute_error関数、mean_absolute_percentage_error関数を使用することにします。

print('RMSE: ',np.sqrt(mean_squared_error(df_test, test_pred)))
print('MAE: ',mean_absolute_error(df_test, test_pred)) 
print('MAPE: ',mean_absolute_percentage_error(df_test, test_pred))

RMSE:  10117.15204491857
MAE:  8592.666666666666
MAPE:  0.020711765602020024

RMSE: 10117

平均的に予測値が実測値から10117人の範囲でずれていることを示します。この値が大きいほどモデルの精度が低いことを示します。
MAE: 8592

平均的に予測値が実測値から8593人の範囲でずれていることを示します。これは誤差の大きさを直接示しており、RMSEと比較すると外れ値の影響を受けにくいです。
MAPE: 0.0207 (約2.07%)

平均的に予測値が実測値から2.07%の範囲でずれていることを示します。MAPEが低いほど予測精度が高いことを示します。

考察

グラフを見ると大きく予測を外しているのにMAPEが2%であることに驚きました。しかし、実測値が毎月40万人前後あることを考えると、8600人の誤差は全体の数%に過ぎないと考えるべきでしょうか。そこで、式(4)、式(5)を用いて平均実測値$\bar{y}$を求めることにしました。

(4)より

$$
\text{MAE} = \frac{1}{n} \sum_{i=1}^{n} \left| y_i - \hat{y}_i \right|=8600\qquad(4)'
$$

(5)より

$$
\text{MAPE} = \frac{100}{n} \sum_{i=1}^{n} \left| \frac{\hat{y}_i - y_i}{y_i} \right|=0.02\qquad(5)'
$$

(4)’を(5)’に代入

$$
\frac{8600}{\bar{y}}=0.02
$$
$$
\bar{y}=\frac{8600}{0.02}=430000
$$

したがって平均実測値$\bar{y}$＝43万人と導くことができました。

実測値のスケールが非常に大きいため、PMSE、MAEが大きくても、割合としては小さいことが確認できました。また、MAEがRMSEよりも小さい値を示しているため、外れ値の影響がある程度抑えられていることがわかります。

今回のSARIMAモデルは、MAPEが2%であることからデータに対して比較的高い予測精度を示していると結論付けます。

今後の展望

都道府県別モデルの構築

全国の総献血者数というスケールが大きすぎるデータでは、地域ごとの特性や違いを反映しにくい可能性があります。そこで、都道府県別にモデルを構築することで、地域ごとの特性を考慮した予測を行います。
説明変数の追加
- 平均降水量や平均気温
  天候は人々の行動に影響を与えるのではないか。
- インフルエンザなどの感染者数
  感染症の流行状況は人々の行動に影響を与えるのではないか。
- 献血ルームの最寄り駅の平均乗降者数
  人の流れが多い場所に位置する献血ルームは、訪れる人の数が多くなるのではないか。
- 「献血」というワードの検索数
  献血に対する関心の高さを示すのではないか。
- 献血イベントの実施回数
  イベントは直接的に献血者数を増加させる要素なのではないか。
新しいモデルの導入
- ニューラルネットワークや深層学習
  単純に興味があることもあるが、ニューラルネットワークや深層学習モデルを活用することで、予測精度を向上させられないか。
- 回帰木
  データの非線形性や交互作用を捉えることに長けている。これにより、SARIMAモデルが捉えにくいパターンを補完することができないか。
- アンサンブル学習
  47都道府県別に作成した複数のモデルを組み合わせたアンサンブル学習を用いることで、予測精度を向上させられないか。

おわりに

今回、授業以外の取り組みとして、初めて時系列分析を一から行い、SARIMAモデルを用いて全国の総献血者数を予測してみました。PDFのOCRから始めて、モデルを構築するまでの過程を通じて、実際のデータを扱いながら問題解決に取り組む良い機会となりました。今後も継続的に学習に励みたいと思います。

GitHubのリポジトリを置いておきます。

宣伝

千葉県学生献血推進協議会
私が所属している千葉県学生献血推進協議会では、自身が献血への知識や理解を深めるための献血セミナーの開催や、年間で５回程の献血啓発イベントなどを行い、献血への理解、興味を持ってもらうような活動を行っています！千葉県内の学生で気になる方がいましたら是非ご参加ください！

参考にしたサイト等

献血者のデータの節
時系列解析 SARIMAモデルを用いた日本のCOVID-19 感染予測実証実験

モデルの構築と予測の節
株式会社セールスアナリティクス｜ARIMA系モデルで予測する方法