Prophetを用いた時系列回帰分析で沼津市の人口推移の将来予測する

やりたいこと

・Pythonの時系列予測ライブラリ「Facebook Prophet」で沼津市の過去数年分の人口動態を元に未来年月の人口を将来予測する
・分析結果を評価⇒チューニング⇒再分析という一通りの流れをやってみる

環境

**　Windows10 64bit**
**　Anaconda3**
　　仮想環境：python(3.7.3)
　　ライブラリ：fbprophet(0.6),pandas(1.1.2),matplotlib(3.0.3)

使用データ

沼津市　「人口と世帯」一覧表

・上記サイトの年度別ファイルから月別の人口を集計したものを入力ファイル**[dataset.csv]**として使用します。
・元々、沼津市のオープンデータには人口統計書があります。
　ふじのくにオープンデータカタログ（沼津市）
　これらは年次統計にあたりますが、時系列分析のデータの粒度としては月次を想定していたので、
　今回は都合よく月次単位のデータが載っていた市公式サイトの住民基本台帳の情報を利用します。
・代償として、PDFからスプレッドシートへのフォーマット整形が必要になりました。面倒でした。
・年月、世帯人員の２列構成です。今回はお試しなので外部変数は無し。
・2013/04～2020/08を学習データとします。

dataset.csv

年月	世帯人員（計）
2013-04	205,790
2013-05	205,739
2013-06	205,572
・・・	・・・
2020-07	193,896
2020-08	193,756

予測分析

###ライブラリのインポート

import pandas as pd
import matplotlib.pyplot as plt
from fbprophet import Prophet

データ取込

Jinkoudata = pd.read_csv("dataset.csv")
Jinkoudata.columns = ['ds','y']

・人口データをデータフレーム[Jinkoudata]に取り込みます。
・列名[年月][世帯人員（計）]をそれぞれ[ds][y]に変更します(prophetの命名規則によるもの)。

モデル生成

model = Prophet(weekly_seasonality=False,daily_seasonality=False)

・今回は月次ベースの周期性にて分析を行うため、週周期・日周期を考慮外とするパラメータを指定。
・その他、インスタンス生成時に下記のようなパラメータを任意に指定できます。
　[]内はパラメータ値。例として指定されている値は規定値。

●growth='linear' [logistic/linear]
　分析モデルを線形（linear）あるいはロジスティック曲線（logistic）で指定。

●changepoints=None [DATE(YYYY-MM-DD)]
　トレンド（長期的な傾向変動）の変化点の任意指定。
　prophetでは変化点を自動抽出して予測に反映しますが、必要に応じて変化点の自動検出を使用する代わりに、
　このパラメータを使用して潜在的なチェンジポイントの場所を手動で指定する意図で用います。

●n_changepoints=25 [任意の定数]
　検出する変化点の候補の数を指定。
　大きければ大きいほど多くの変化点が検出されやすくなります。
　これを大きくし過ぎると変化点が過剰に検出され、予測精度を落とす場合があります。

●yearly_seasonality='auto' [true/false/auto/任意の定数]
　周期性（年）を指定。任意の数字を指定することで考慮度合いを調整可能（指定がない場合は10が規定値）。

●weekly_seasonality='auto' [true/false/auto/任意の定数]
　周期性（週）を指定。任意の数字を指定することで考慮度合いを調整可能（指定がない場合は10が規定値）。

●daily_seasonality='auto' [true/false/auto/任意の定数]
　周期性（日）を指定。任意の数字を指定することで考慮度合いを調整可能（指定がない場合は10が規定値）。

●holidays=None [DataFrame]
イベント日付のデータフレームを指定。
　祝祭日やイベントのある日付を分析に組み込み、予測に反映する。

●changepoint_prior_scale=0.05 [任意の定数]
　ラプラス分布の分散を指定。この値が大きければ大きいほど変化点が検出されやすくなります。
　しかしそれにより、柔軟性が高くなりすぎる場合があります。n_changepointsと同様。

●changepoint_range=0.8 [任意の定数]
　変化点検出を行うデータの割合を指定。
　規定値は0.8なので、学習データの先頭から８割分を対象に変化点を抽出しますが、
　例えばデータの終盤でトレンドに大きく動きがある場合に「1.0」と設定することで、
　終盤の変動も含めて全範囲のデータを使って変化点を検出するようになります。

　　他にも色々あるようですが、今のところの知識量で考慮できそうなのが上記くらいなので割愛します。

学習

model.fit(Jinkoudata)

fit()
・fitで学習データを指定することで学習します。

データフレーム生成

future = model.make_future_dataframe(periods=16,freq='M')

make_future_dataframe()
・作成したモデルに学習データとしてfitしたデータフレームに加える形で、未来の時系列行を作成します。
・ここでは、fitしたデータフレーム（Jinkoudata：2013/04～2020/08）に加えて、
　その先の16ヶ月分（2020/09～2021/12）のデータフレームを作成しました。
・結果として、[future]は2013/04～2021/12という期間を持つデータとなります。

future

	ds
0	2013-04-01
1	2013-05-01
2	2013-06-01
・・・	・・・
103	2021-10-31
104	2021-11-30

予測

predict = model.predict(future)

・先ほど作成した時系列のデータフレームを引数に指定してpredictします。
・予測対象となる時系列はfutureに格納されている時系列である2013/04～2021/12となります。

予測結果 [predict]

	ds	yhat	yhat_lower	yhat_upper
0	2013-04-01	205798.225901	205762.014035	205838.079234
1	2013-05-01	205717.024368	205679.523657	205759.642408
2	2013-06-01	205600.30351	205560.587079	205640.792342
3	2013-07-01	205530.899822	205488.087052	205570.056443
4	2013-08-01	205335.468464	205296.411528	205378.093757

予測（predict）してデータフレームに格納される情報は他にもありますが、
数値の捉え方が難しいためここでは予測値関連の項目を確認します。
　・yhat：予測値
　・yhat_lower：予測下限値
　・yhat_upper：予測上限値
上記の列に、「予測値」「予測下限値」「予測上限値」が格納されます。
「予測値」は可能性の高い予測結果、「予測下限値」「予測上限値」は予測値の振れ幅としての下限・上限です。

分析

model.plot(predict)

・予測の描画は[plot]で行います。
・黒い点が学習データ、青い線が予測値（yhat）、青い帯が予測の振れ幅（yhat_lower：yhat_upper）です。

【ここから読み取れること】
　■予測結果として、2021年についてもこれまでと同様のペース（およそ1000人以下/年）での減少傾向
　■2021年中に193000人台を割る
　■予測値と実績は離れていないことから予測の精度は高め（というよりデータが素直）
　■およそ半年先までは比較的予測値の振れ幅の少ない予測が可能
　　（学習データは2020/08まで、2021/03以降から振れ幅が急に拡大している）
　

model.plot_components(predict)

・傾向変動、季節変動のような因子毎の描画は[plot_components]で行います。
・[trend]項がトレンド（長期的な傾向変動）。年単位の推移としてどういう変動であるかを描画しています。
・[yearly]項が季節変動。１年の中での変動を描画しています。
・トレンドは、短期間的な増減はあまり表さず、季節変動を繰り返していく１年毎の規模増減を観測すます。

【ここから読み取れること】
　■ここ８年程、綺麗な線形状に減少中。学習開始年次（2013年）以降、年単位で盛り返す兆候無し
　■年内の変動として、基本的には減少が続くが、１０・１１月に僅かに増加に転じる
　　（季節変動は波形で分析・表現されるため、今回のような月次データを用いている場合に
　　　データの日付間隔が大きい分、波形が大きくなることが想定されるので注意が必要）
　　（実データ上は増加はほぼないが、１０・１１月は他より減り幅が顕著に縮小してきていることと、
　　　近年微増を見せたていることから予測として増加に転じる箇所が生まれたと思われる多分）
　■僅かながら、減少傾向は緩やかになりつつある
　　（例年減少が激しいのは３月だが、減り幅自体も縮小傾向にある。
　　　例えば３月の減少自体は就職・進学のイベントによるもの、
　　　減り幅の縮小は少子化によるものという推測もできるが、
　　　詳細を詰めるには転入出のような社会動態、出生・死亡者数のような
　　　自然動態といったオープンデータが必要）

評価

作成したモデルがどの程度の予測精度を有しているのか、評価が必要になります。
今回は交差検証を用いて評価を行っていきます。

交差検証

from fbprophet.diagnostics import cross_validation
df_cv = cross_validation(model, initial='1095 days', period='182 days', horizon = '365 days')

cross_validation()

　・initial :予測を開始する時期
　・period :データをずらす期間
　・horizon :予測する期間
上記宣言例である場合、1095日目を始点として（initial）、
データの1日～1094日分を学習データ期間（指定がなければhorizonの３倍）、
1095日目から365日分（horizon）を評価期間として予測します。
予測の分析が終わったら予測期間の始点を182日後ろにずらし（period）、
1277日目から365日分を次の評価期間として予測します。
これをデータの終わりまで繰り返し、予測結果を格納します。
この予測結果は、後述の方法で確認しやすい尺度に置き換えて評価結果とします。

交差検証結果 [df_cv]

	ds	yhat	yhat_lower	yhat_upper	y	cutoff
0	2016-09-01	199454.428910	199389.987532	199514.165966	199386	2016-08-05
1	2016-10-01	199437.182903	199145.883231	199744.554272	199236	2016-08-05
2	2016-11-01	199274.971653	198645.490967	199906.989114	199110	2016-08-05
3	2016-12-01	199073.586270	198066.379455	200065.689963	199006	2016-08-05
4	2017-01-01	198809.990576	197391.138491	200210.975003	198837	2016-08-05

　

評価尺度計算

from fbprophet.diagnostics import performance_metrics
df_p = performance_metrics(df_cv)

performance_metrics()

交差検証結果[df_cv]に対して、評価尺度を計算し、結果を格納します。

評価尺度計算結果 [df_p]

	horizon	mse	rmse	mae	mape	mdape	coverage
0	57 days	16078.849437	126.802403	109.436237	0.000554	0.000432	0.1250
1	58 days	20194.085707	142.105896	129.175240	0.000655	0.000627	0.1250
2	59 days	18026.561842	134.263032	113.963954	0.000578	0.000545	0.2500
3	60 days	20932.092431	144.679274	129.291089	0.000657	0.000687	0.2500
4	87 days	33783.385041	183.802571	161.495031	0.000819	0.000824	0.3750

評価尺度は、結果を評価指標別に計算したもので、種類としては以下の通り。

●MSE [平均平⽅誤差(mean squared error)]
　実際の値と予測値の差となる絶対値の二乗を平均した値。

●RMSE [平均平⽅誤差の平⽅根(root mean squared error)]
　実際の値と予測値の差の二乗の平均値の平方根。

●MAE [平均絶対誤差(mean absolute error)]
　実際の値と予測値の差の絶対値の平均。

●MAPE [平均絶対パーセント誤差(mean absolute percent error)]
　実際の値と予測値の差の実際の値に対する割合の平均。

MSE、RMSE、MAEはそれぞれ算出は異なりますがいずれも予測のずれを平均的に実測値ベース（今回は人口）で表現します。
MAPEは予測のずれを割合ベース（％）で表現します。
今回は扱う数字も大きいため、感覚的にわかりやすいMAPE（平均絶対パーセント誤差）にて評価結果を確認したいと思います。

MAPE（平均絶対パーセント誤差）による評価尺度可視化

from fbprophet.plot import plot_cross_validation_metric
fig = plot_cross_validation_metric(df_cv, metric='mape')

plot_cross_validation_metric()

前項で計算した評価尺度を可視化できます。
交差検証にて指定したHorizon（評価期間）をＸ軸に、
実行された評価回数分のMAPEの値が点で描画され、
その平均が青線で表示されている形です。

【ここから読み取れること】
　■全体的な予測の精度は低くない
　　⇒最大の外れ値を見てもmapeは0.4%程。
　　（＝20万人規模の自治体であることから、概算で最大800人相当のずれ）
　　　平均を辿ると0.05～0.15%の間で誤差が推移しており、100～300人程のずれです。
　　　あくまで予測なので正解とする指標はありませんが、
　　　実測データ上は月次で100人単位の動きがある期間も多いため、
　　　誤差としても１集計期間で埋まる程度の差であることから
　　　大幅な誤差ではないという判断ができます。
　■予測開始から期間が経過するほど精度が落ちる
　　⇒当然ながら月が経つほど精度が落ちる（mapeの分散が大きくなる）傾向があります。

チューニング

モデル精度向上の余地

作成したモデルが最善であるとは限らないため、チューニングを施して作成したモデルを用いて予測することで精度の向上を図る場合があります。
ここでは「チューニングの余地があるか」といった判断のプロセスも併せて書いていきます。

チューニングの観点

チューニングを施すか否かといった有無を判断する場合、学習に用いたデータの性質を元に考える必要があります。
今回は下記のような観点で考えていきます。

　●加法的季節性・乗法的季節性
　　実データの変動は加法的（目的変数が大きさに関わらず変動の幅・大きさは一定）なものであるか、乗法的（目的変数が大きければ変動の幅も大きくなる）なものであるか

　●変化点の抽出数
　　実データの変化点の個数は適切か

　●変化点の抽出範囲
　　実データの変化点の抽出範囲は適切か

今回のモデルへのチューニングの判断

　●加法的季節性・乗法的季節性
　　今回の目的変数は人口データです。
　　人口増減の１要因と考えられる「死亡」「出生」は高齢層や若年層といった
　　人口の層の母数により変動の大きさも変化すると考えられるため、
　　今回のデータは僅かながらに乗法的な一面を持っていると判断します。
　　（予測結果のplotからも、年々全体数が少なくなっていくもののそれに伴って
　　　ほんの少しですが年単位の減少が小さくなっていることが確認できます）
　　（他には「転入」「転出」も同様に人口の母数の影響を受けますが、
　　　これらの場合は同時に転入出の受け口となる企業・学校の数や規模の影響も色濃く受けるため、
　　　人口データのみでの断定は難しかったです）
　　Prophetのモデル作成はパラメータ[seasonality_mode]で「加法的季節性」または
　　「乗法的季節性」を指定できますが、省略時はデフォルトで「加法的季節性」でモデル作成します。
　　チューニングとして「乗法的季節性」でのモデル作成を試してみる余地があります。
　　⇒チューニング対象

　●変化点の抽出数
　　Prophetでは学習データの中から変化点を自動で抽出して予測に用いており、
　　抽出する変化点の個数もデフォルトで25となっています。
　　変化点の抽出数はパラメータ[n_changepoints]でモデル作成時に指定が可能です。
　　この数が少なすぎると「この時期は変化が起きる」といった
　　時系列上の変化点の認識が不足するため予測の精度が落ちますが、
　　この数が多すぎても僅かな変動を変化点として認識してしまい、
　　予測の柔軟性を損ないやはり精度が落ちる想定です。
　　つまり数を指定するのであればモデルに学習させたい
　　肝要な変化点を押さえる必要がありますが、
　　今回のような月次毎のデータではイベントが不明瞭になりがちで
　　「何個が適正か」といった判断が難しいため、チューニングの余地は薄いと考えます。
　　⇒チューニング対象外

　●変化点の抽出範囲
　　Prophetでは変化点の抽出範囲として、
　　学習データの始点から80%を対象期間として抽出します。
　　例えば、学習データの終盤に急激なトレンド上昇があった場合、
　　デフォルトのモデルではその変化を考慮した予測が行えませんので、
　　モデル作成時にパラメータ[changepoint_range]にて
　　'1.0'を指定するなどして抽出範囲を調整します。
　　今回の場合はそういった直近での大きな変動といったものは
　　特に見受けられないため、変化点の抽出範囲という観点でのチューニングの余地はなさそうです。
　　⇒チューニング対象外

　

チューニング：乗法的季節性でのモデル作成

上項にて判断した通り、チューニングとして「乗法的季節性」でのモデル作成を行っていきます。

# モデル生成（乗法的季節性）
# 生成時パラメータを追加 [seasonality_mode = "multiplicative"]
model_multi = Prophet(weekly_seasonality=False,daily_seasonality=False,seasonality_mode = "multiplicative")
# 学習
model_multi.fit(Jinkoudata)
# 予測
predict_multi = model_multi.predict(future)
# 予測結果描画
fig = model_multi.plot(predict_multi)

乗法的モデルの予測結果の描画は下の通りです。


傾向としては加法的モデルと同様に見えますね。
こちらについても交差検証をしていきます。

# 交差検証
df_cv_multi = cross_validation(model_multi, initial='1095 days', period='182 days', horizon = '365 days')
# 評価尺度計算
df_p_multi = performance_metrics(df_cv_multi)
# MAPE（平均絶対パーセント誤差）による評価尺度可視化
fig_multi = plot_cross_validation_metric(df_cv_multi, metric='mape')

評価尺度MAPEは以下の通りでした。

僅かに、本当に極僅かにですが、加法的モデルより分散が少なく、
予測とのずれが少ないモデルとなりました。

実測データと２種の予測結果を同じ表に描画して比較してみます。

# 表サイズ設定
plt.figure(figsize=(12,6))
# 軸見出し設定
plt.xlabel('Year month')
plt.ylabel('Jinkou')
# 描画範囲指定
plt.xlim(82, 105)
plt.ylim(190000, 196000)
# 目盛りの表示変更
plt.xticks([82,88,94,100,105], ["2020-01","2020-07","2021-01","2021-07","2021-12"])
# 各数値の描画
plt.plot(Jinkoudata.y,label = 'Jinkoudata')
plt.plot(predict.yhat,label = 'Predict_add')
plt.plot(predict_multi.yhat,label = 'Predict_multi')
plt.legend()

描画結果は下の通りでした。

[Jinkoudata]が学習に用いた実測データ、
[Predict_add]が加法的季節性モデルでの予測結果、
[Predict_multi]が乗法的季節性モデルでの予測結果です。
……これでいいのか、という疑問はありますが、[Predict_add]と[Predict_multi]の線は重なっておりほぼ同一といっていい結果になりました。
元々、実測データの乗法的な性質自体は僅かなものだったので、予測結果としてもさほどの差として現れないのかもしれません。

まとめ

モデル作成や予測自体は思った以上に簡単でした。
ただし活用にはやはりデータの性質や相関性についての理解・知識が必要になりますね。
将来予測のモデルとして信頼していい予測期間というのは思ったより長くなかったので、（今回の例に限らなさそうですが）単位期間毎にデータから傾向を読み取って都度チューニングする、というのは運用として前提になりそうかなと思いました。

練習ということで人口の将来予測をしてみましたが、元のデータからして減少推移は明らかだったので、ここまでやっておいてなんですがあまり予測の価値はなかったと思います。テーマ選びはやはり大事ですね。
深堀りするなら外部変数を使って精度を上げる等が考えられますが、その相関分析を別途試してみるのもいいかもしれません。
沼津市は好きな街なので別の分野も勉強がてら予測してみるのもいいかもと思いました。