2023年 日本出生数予測と実際の結果について
はじめに
貴重なお時間ありがとうございます。
Aidemy Premiumにて「データ分析コース」を６ヶ月受講いたしました。
Aidemy Premiumを受講したのは、データ分析、機械学習、Pythonをしっかりと学びそれを活かして仕事をしていきたいと思ったことが、きっかけでした。
AI、機械学習、ビックデータ、データ活用とたくさんのデータに関するワードを耳にするようになりました。
そしてデータサイエンティストという新たな仕事に就く人が不足していること、大学にデータサイエンスに関する学部や学科、専攻が増えてきていることを知りました。
その中で興味を持ち、数学が得意なこともありデータサイエンスについて学習することにしました。
ただ、実際に調べてみると難しい用語やコードが多くなかなか理解できないケースが多くあるかと思います。
記事を書く１番の目的は知識を定着させることですが、
少しでも皆様と知識や経験をシェアできるようにレベルをあげていきます。
実際にPythonで実装した内容で記載しております。
わかりにく表現もあるかと思いますが、よろしくお願いいたします。
実行環境
・GoogleColabolatry
・Python 3.11.7
データセット
Kaggleにある日本の出生人口統計よりダウンロード。
※ 1899 年から 2022 年までの日本の出生関連統計の集合データ。第二次世界大戦中に記録が失われたため、1944 年から 1946 年までの一部のデータが欠落しています。
日本の出生数について
2023年の出生数が2024年2月に出生数75万8631人と厚生労働省が速報値として発表しました。
そこで、2022年までの総人口や出生数のデータをもとに2023年の出生数を予測し、実際の出生数と比較することにしました。
まず、最初に時系列分析をかけ、予測値を求めました。
次に単回帰分析や重回帰分析を行うために、ヒートマップを表示して、出生数と相関の高いデータを探し、１つまたは２つの項目に絞って分析にかけました。
その数値をもとに分析や実験を繰り返し行いました。
データ前処理について
ライブラリ
import numpy as np
import pandas as pd
import seaborn as sns
import missingno as msno
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib.ticker import ScalarFormatter
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
import matplotlib.pyplot as plt
import warnings
import itertools
from sklearn.linear_model import LinearRegression
from sklearn.datasets import make_regression
from sklearn.model_selection import train_test_split処理
124行50列のデータになります。1944年から1946年の一部データが欠損値となっています。
df = pd.read_csv("/kaggle/input/japan-birth-statistics/japan_birth.csv", index_col=0)
formatter = mpl.ticker.StrMethodFormatter('{x:,.0f}')print(df.columns.values)
type(df.columns.values)各カラムを出力してみます。
['year' 'birth_total' 'birth_male' 'birth_female' 'birth_rate'
'birth_gender_ratio' 'total_fertility_rate' 'population_total'
'population_male' 'population_female' 'infant_death_total'
'infant_death_male' 'infant_death_female' 'infant_death_unknown_gender'
'infant_death_rate' 'infant_death_gender_ratio'
'infant_deaths_in_total_deaths' 'stillbirth_total' 'stillbirth_male'
'stillbirth_female' 'stillbirth_unknown_gender' 'stillbirth_rate'
'stillbirth_gender_ratio' 'firstborn' 'secondborn' 'thirdborn'
'forthborn' 'fifthborn_and_above' 'weeks_under_28' 'weeks_28-31'
'weeks_32-36' 'weeks_37-41' 'weeks_over_42' 'mother_age_avg'
'mother_age_firstborn' 'mother_age_secondborn' 'mother_age_thirdborn'
'mother_age_under_19' 'mother_age_20-24' 'mother_age_25-29'
'mother_age_30-34' 'mother_age_35-39' 'mother_age_40-44'
'mother_age_over_45' 'father_age_avg' 'father_age_firstborn'
'father_age_secondborn' 'father_age_thirdborn' 'legitimate_child'
'illegitimate_child']出力結果を最初の5行だけ載せます。列は最初の5列(年、出生数、男の子の出生数、女の子の出生数、出生率)になります。
year birth_total birth_male birth_female birth_rate \
0 1899 1386981.0 713442.0 673539.0 32.0
1 1900 1420534.0 727916.0 692618.0 32.4
2 1901 1501591.0 769494.0 732097.0 33.9
3 1902 1510835.0 773296.0 737539.0 33.6
4 1903 1489816.0 763806.0 726010.0 32.7 これらの中から必要なデータを用いて各分析を行っていきたいと思います。
df2=df
ax = df2.plot(x="year", y=["birth_total"], figsize=(12,5))
ax.yaxis.set_major_formatter(formatter)
plt.show()横軸は年度(year)、縦軸は出生数(birth_total)になります。
分析手法
・時系列分析
まずはグラフにし確認していきます。
・原系列
なにもされていない時系列データそのものを 原系列と呼びます。時系列分析において原系列そのものを扱うことはあまりありません。実際には時系列データを加工し新しい系列にして、それを分析してモデルを構築していきます。
・時系列データのパターン
時系列データには トレンド、 周期変動、 不規則変動の3つのパターンがあります。
1. トレンド
データの長期的な傾向を意味します。時間の経過とともにデータの値が上昇または下降している時系列データは「トレンドがある」といい、値が上昇している場合には正のトレンド、下降している場合には負のトレンドがあると言います。
2. 周期変動
周期変動がある場合、時間の経過に伴ってデータの値が上昇と下降を繰り返します。特に1年間での周期変動を季節変動といいます。
3. 不規則変動
時間の経過と関係なくデータの値が変動することをいいます。
年度ごとの出生数をもとに時系列分析にかけました。
データ処理については以下の通りです。
出生数(birth_total)の列と年度(index)に1899年から2022年をもとに時系列
分析します。
birth_total = df["birth_total"]
index = pd.date_range('1899-12-31', '2022-12-31', freq='Y')
birth_total.index = index 時系列解析で欠かせないものの一つに定常性があります。
定常性とは、ある時系列データの値ｙが、どの時点ｘにおいても平均、
分散、自己共分散が一定であるという性質のことです。
これは弱定常性とも言います。
簡単にいうと、各時点での値がある値を中心にばらついているということ
です。
SARIMAモデルとは、自己回帰モデル(ARモデル)と移動平均モデル(MAモデ
ル)と和分モデル(I)を組み合わせた自己回帰和分移動平均モデル(ARIMAモ
デル)に季節的な周期変動を取り入れたモデルです。
SARIMAモデルの理解には上記のモデルの理解が必要不可欠なので、
まずはそれらのモデルについて説明します。
1. ARモデル
ARモデルは、自身の過去の値で回帰を行うことで現在の値を予測するとい
うモデルになります。
2. MAモデル
MAモデルは、過去の誤差の積み重ねにより現在の値を求めるというモデ
ルになります。
3. ARMAモデル
ARMAモデルはARモデルとMAモデルを組み合わせたモデルになります。
4. ARIMAモデル
時系列データによっては、その系列が定常性を満たさず、階差をとった系
列が定常性を満たすものもあります。そのような時系列データについては
ARIMAモデルを使用することができます。
p は 自己相関度といい、 モデルが直前p個の値を用いて予測されるのかを
表す。
d は 誘導といい、 時系列データを定常にするためにd次の階差が必要だっ
たこと を表す。
q は 移動平均といい、 モデルが直前q個の値に影響を受けること を表す。
SARIMAモデル とはARIMAモデルをさらに季節周期を持つ時系列データに
も拡張できるようにしたモデルです。SARIMAモデルは(p, d, q)のパラメー
ターに加えて(sp, sd, sq, s)というパラメーターも持ちます。
sp,sd,sq はそれぞれ 季節性自己相関、 季節性導出、 季節性移動平均とい
い、現在のデータはひとつ以上の季節期間を経た過去のデータに影響され
ます。
例えば12ヶ月周期の季節変動を持つデータの場合、 s のパラメーターは周
期を表すため、s=12 となります。
sq=1 ならちょうど12ヶ月(1周期前)のデータの影響をモデルが受けること
を表します。
時系列データが定常性を満たしているかどうかを確認するためには、視覚的な分析と統計的な検定の二つの方法が一般的に用いられます。
今回はADF検定を用いました。
ディッキー-フラー検定（ADF検定）: ディッキー-フラー検定は、単位根検定の一つで、時系列データが非定常である（単位根を持つ）という帰無仮説を検定します。p値が一定の有意水準（例えば0.05）以下であれば、帰無仮説を棄却し、データが定常であると結論づけます。
def adf_test(df2):
adf_df = pd.DataFrame(
[
adfuller(df2)[1]
],
columns=['P値']
)
adf_df['P値'] = adf_df['P値'].round(decimals=3).astype(str)
print(adf_df)
index = pd.date_range('1899-12-31', '2022-12-31', freq='YE')
# 原系列のp値を算出
adf_test(df2=index)今回の結果はp=0.954でした。そのため、非定常であると結論づけました。
そこで、原系列、移動平均、原系列から移動平均を引いたもの、階差とグラフにしてみました。
df2 = df["birth_total"]
plt.subplot(6, 1, 1)
plt.xlabel("year")
plt.ylabel("birth_total")
plt.plot(df2)
# 移動平均を求める
df2_move = df2.rolling(window=5).mean()
# 移動平均のグラフ
plt.subplot(8, 1, 3)
plt.xlabel("year")
plt.ylabel("birth_total")
plt.plot(df2_move)
# 原系列-移動平均グラフ
plt.subplot(8, 1, 5)
plt.xlabel("year")
plt.ylabel("birth_total")
df2_sabun = df2-df2_move
plt.plot(df2_sabun)
plt.show()
# 階差
df2_kaisa = df2.diff()
plt.subplot(8, 1, 7)
plt.xlabel("year")
plt.ylabel("birth_total")
plt.plot(df2_kaisa)
plt.show()すると、徐々に定常過程近づいていることがわかります。
よって、SARIMAモデルにて2023年の出生数の予測を行いました。
今回は年度ごとのデータであるため、季節変動が見られないため、一旦周
期を12として分析しました。
birth_total = df["birth_total"]
index = pd.date_range('1899-12-31', '2022-12-31', freq='Y')
birth_total.index = index
def selectparameter(DATA, s):
p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], s) for x in list(itertools.product(p, d, q))]
parameters = []
BICs = np.array([])
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod = sm.tsa.statespace.SARIMAX(DATA,
order=param,
seasonal_order=param_seasonal)
results = mod.fit()
parameters.append([param, param_seasonal, results.bic])
BICs = np.append(BICs, results.bic)
except:
continue
return parameters[np.argmin(BICs)]
best_params = selectparameter(birth_total, 12)
SARIMA_birth_total = sm.tsa.statespace.SARIMAX(birth_total, order=best_params[0],
seasonal_order=best_params[1],
enforce_stationarity=False, enforce_invertibility=False).fit()
pred = SARIMA_birth_total.predict('2023-01-31', '2024-12-31')
print(pred)
plt.plot(birth_total)
plt.plot(pred, "r")
plt.show()結果として2023年の出生数の予測値は73万6403人でした。
・単回帰分析
まずはデータのヒートマップを調べてみました。
すると、出生数(birth_total)に対して、出生後にわかる数値ばかりであった
ため、今回は総人口から予測してみることにしました。
key_columns = ['year', 'birth_total', 'birth_rate',
'total_fertility_rate', 'population_total',
'infant_death_total', 'infant_death_rate',
'infant_deaths_in_total_deaths',
'stillbirth_total', 'stillbirth_rate',
'weeks_under_28', 'weeks_28-31',
'weeks_32-36', 'weeks_37-41', 'weeks_over_42', 'mother_age_avg',
'mother_age_under_19', 'mother_age_20-24', 'mother_age_25-29',
'mother_age_30-34', 'mother_age_35-39', 'mother_age_40-44',
'mother_age_over_45', 'father_age_avg', 'legitimate_child',
'illegitimate_child']
key_df = df[key_columns]
plt.figure(figsize=(15,10))
sns.heatmap(key_df.corr(),annot=True,cmap='Reds')
plt.show()
2022年の人口は調べてみたところ、10月の人口を用いていました。
2023年の人口が分析時に発表がされておらず、9月の人口までしか発表
されていなかったため、今回は9月の人口をもとに予測してみました。
df.dropna(inplace=True)
X = df[["population_total"]]
y = df[["birth_total"]]
train_X, test_X, train_y, test_y = train_test_split(X, y, random_state=42)
a = [[121270000]]
a = np.array(a)
# モデルの構築
model = LinearRegression()
# モデルの学習
model.fit(train_X,train_y)
# test_X, test_yに対する決定係数を出力してください
print(model.score(test_X, test_y))
print(model.predict(a)) 2023年の人口1億2127万人に対して、単回帰分析を行ったところ、出生
数の予測値として86万9131人となりました。
・重回帰分析
同じようにヒートマップを参考にし、男性と女性の人口から重回帰分析を
することにしました。
単回帰分析の時と同様に、2023年の人口が発表されていないため、9月
時の男性と女性の人口を用いました。
男性が5893万7000人、女性が6233万3000人となります。
X = df[["population_male","population_female"]]
y = df[["birth_total"]]
model = LinearRegression()
model.fit(X, y)
a = [[58937000,62333000]]
a = np.array(a)
train_X, test_X, train_y, test_y = train_test_split(X, y, random_state=42)
# モデルの構築
model = LinearRegression()
# モデルの学習
model.fit(train_X, train_y)
model.score(test_X, test_y)
print(model.predict(a))その結果、2023年の出生数の予測値は74万4214人でした。
分析方法の比較
以上の結果それぞれの手法による予測値は以下の通りです。
時系列分析 73万6403人
単回帰分析 86万9131人
重回帰分析 74万4214人
そしてニュースで発表された速報値は75万8631人でした。
単純に見ると、重回帰分析が一番結果に近い予測を出したことになります。
しかし、一般的に回帰分析よりも時系列分析の方が、より正確な予測値を出してくれることがわかっています。
では、なぜこのような結果がでたのでしょうか。
これは、月ごとのデータではないため、今回時系列分析の周期を12としたことが理由であると考えられます。
そのため、周期をいくつにすればより速報値に近い値がでるのか調べてみました。
s= 2 73万7464人
s= 3 74万4944人
s= 4 74万3724人
s= 5 75万5261人
s= 6 74万869人
s= 7 74万9267人
s= 8 77万1755人
s= 9 74万3450人
s=10 75万8275人
s=11 74万8391人
s=12 73万6403人
s=13 77万5725人
s=14 78万3635人
s=15 75万6973人
s=16 75万42人
s=17 75万3087人
s=18 72万5006人
s=19 67万4318人
s=20 76万6943人
s=20まで調べたが、その中で一番速報値に近い値となったのはs=10であった。s=21以降は調べなかったが、sの値が5の倍数の時に速報値にやや近い値になる傾向があるので、より近い値が求められるとしたら5の倍数ではないかと考えます。
最後に
データ分析初心者でしたが、分析結果を今回まとめました。もし、訂正すべき点や助言などいただければ幸いです。
つたない文章でしたが、最後までお読みいただきありがとうございました。
この記事が気に入ったらサポートをしてみませんか？