日本の超過死亡?

厚労省人口動態調査から,人口動態統計(確定数),人口動態統計月報(概数),人口動態統計速報がリンクされている。確定数は翌年9月ごろにならないと出ない。概数はほぼ5ヶ月後,速報はほぼ2ヶ月後に出る。

データはe-Statの人口動態調査で公開されている。ここから集めた確定数・概数・速報データを japandeaths.csv というファイルに収めた。確定数は上記ページから確定数→死亡→年次→2020年とたどって,5-4「死亡月別にみた年次別死亡数及び死亡率(人口千対)」を選ぶ。このCSVファイルは整然としていないのでデータの整然化の例のようにして整然化する。

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import datetime

df = pd.read_csv("https://okumuralab.org/~okumura/python/data/japandeaths.csv")
t = [datetime.datetime(int(x["年"]), int(x["月"]), 1) for i, x in df.iterrows()]

def days(year, month=None):
    year = int(year)
    leap = (year % 4 == 0) and (year % 100 != 0) or (year % 400 == 0)
    if month is None:
        return 365 + leap
    else:
        month = int(month)
        m = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
        return m[month - 1] + ((month == 2) and leap)

perday1 = np.array([r[2] / days(r[0], r[1]) for i, r in df.iterrows()]) # 確定数
perday2 = np.array([r[3] / days(r[0], r[1]) for i, r in df.iterrows()]) # 概数
perday3 = np.array([r[4] / days(r[0], r[1]) for i, r in df.iterrows()]) # 速報値

plt.plot(t, perday1, "o-")
plt.plot(t, perday2, "s-")
plt.plot(t, perday3, "*-")
plt.ylabel("1日あたり死亡数")
plt.legend(["確定数", "概数", "速報値"])
Deaths in Japan
df["死亡数"] = [n3 if np.isnan(n1) and np.isnan(n2) else n2 if np.isnan(n1) else n1
                for n1, n2, n3 in zip(df["確定数"], df["概数"], df["速報値"])]
df["日数"] = [days(year, month) for year, month in zip(df["年"], df["月"])]

for y in sorted(set(df["年"])):
    df1 = df[df["年"] == y]
    plt.plot(df1["月"], df1["死亡数"] / df1["日数"],
             alpha=0.5, marker=f"${y % 10}$", label=y)
plt.xlabel("月")
plt.ylabel("1日あたり死亡数")
plt.legend(loc=(0.38, 0.60), labelspacing=0.1)
Deaths in Japan
for m in range(1, 13):
    df1 = df[df["月"] == m]
    plt.plot(df1["年"], df1["死亡数"] / df1["日数"], marker=f"${m}$")
plt.ylabel("1日あたり死亡数")
plt.xticks(sorted(set(df["年"] // 2 * 2).intersection(set(df["年"]))))
Deaths in Japan
df1 = df.groupby("年")[["死亡数", "日数"]].sum()
df2 = df1[df1.index < 2022]
y = df2.index
x = df2["死亡数"] / df2["日数"]

plt.plot(y, x, "o-")
plt.xlabel("年")
plt.ylabel("1日あたり死亡数")
Deaths in Japan

非常に大雑把に2020〜2021年の超過死亡数を求めてみよう。2011年は東日本大地震で死亡数が増えたと考えられるので,ベースラインを2012〜2019年に設定して,最小2乗法で直線をフィットし,そこからの外れを調べよう。

u = (2012 <= y) & (y <= 2019)
slope, intercept = np.polyfit(y[u], x[u], 1)
plt.clf()
plt.plot(y, x - (slope * y + intercept), "o-")
plt.axhline(linewidth=0.5, color="black")
plt.xlabel("年")
plt.ylabel("1日あたり超過死亡数")

厚労省オープンデータによれば,COVID-19による2020年12月31日までの累積死亡数は3459人,2021年12月31日までの累積死亡数は18385人である。これをオレンジ色で描き足す。

plt.plot([2020, 2021], [3459 / 366, (18385 - 3459) / 365], "o-")
Excess Deaths in Japan

このように,2020年の死亡数は予想よりかなり少ない(-68.6人/日)。その反動か,2021年には予想よりかなり多く(+69.3人/日),報告された死亡数(40.9人/日)よりも多い(1.7倍,+28.4人/日)。一部で騒がれた「報告数の6倍」はどう考えてもありえないし,増分(+28.4人/日)は2012〜2019年のゆらぎと同程度なので,何ともいえない。

ちなみに月ごとの1日あたりのCOVID-19による死亡数は次の通りである。

1日あたりCOVID-19死亡数
2020-062.7
2020-071.2
2020-089.2
2020-099.2
2020-106.3
2020-1112.4
2020-1242.6
2021-0172.9
2021-0277.3
2021-0341.1
2021-0435.6
2021-0590.9
2021-0657.7
2021-0713.2
2021-0827.4
2021-0953.6
2021-1020.0
2021-113.1
2021-121.1
2022-0112.9
2022-02172.9
2022-03144.0
2022-0449.1
2022-0533.8
2022-0619.0
2022-0742.1
2022-08235.3
df = pd.read_csv("https://covid19.mhlw.go.jp/public/opendata/deaths_cumulative_daily.csv",
                 parse_dates=['Date'])
a = np.arange(np.datetime64("2020-06"), np.datetime64("2022-10"))
da = pd.to_datetime(a)
da1 = da - np.timedelta64(1,"D")
da1s = pd.Series(da1, name='Date')
df2 = pd.merge(da1s, df).diff()
x = df2.ALL / (df2.Date / np.timedelta64(1,"D"))
for i in range(1, len(x)):
    print(f"{a[i-1]},{x[i]:3.1f}")

Last modified: