超過死亡率の集計

コロナ禍で超過死亡が話題になる事がありましたがどうやって計算しているのかよく知らなかったので少しひも解いてみました。参照した文献はこちら。使用したのはexcessmortパッケージです。

このパッケージにはCDCが公開する、州別で、週ごとの死亡者数の見積もり、に加えて人口の見積もり（Censusという人口統計からの変化を計算した見積もり）の情報がcdc_state_countsというデータフレームで含まれます。これを利用して、新型コロナが流行し始めた2020年以降データと、死亡の多かった2017-18のインフルエンザシーズンを除外したデータを基に、それらの死亡の影響がない場合の期待される死亡のトレンド・季節変化を計算します。そのうえで、そのモデルからの逸脱状況をプロットします。（その逸脱分が新型コロナ（と2017-18はインフルエンザ）の影響と解釈したいという仮説で集計します）

メリーランド州を選択したのは、以前住んでいたからです。

library(excessmort)
library(dplyr)
library(lubridate)
library(ggplot2)



# A list with six ggplot objects: ‘population‘ is a time series plot of the population. ‘seasonal‘ is a
# plot showing the estimated seasonal effect. ‘trend‘ is a time series plot showing the estimated trend.
# ‘weekday‘ is a plot of the estimated weekday effects if they were estimated. ‘expected‘is a time
# series plot of the expected values. ‘residual‘ is a time series plot of the residuals.

flu_season <- seq(make_date(2017, 12, 16), make_date(2018, 1, 16), by = "day")
exclude_dates <- c(flu_season, seq(make_date(2020, 1, 1), today(), by = "day"))
res <- cdc_state_counts %>%
  filter(state == "Maryland") %>%
  compute_expected(exclude = exclude_dates,
                   keep.components = TRUE)

p <- expected_diagnostic(expected = res, alpha = 0.50)
p$population
p$seasonal
p$trend
p$expected
p$residual

このRのコードは、統計解析のために使用されるパッケージや関数を活用しています。以下に、コードの動作をプログラム初心者向けに説明します。

ライブラリの読み込み:
- excessmort, dplyr, lubridate, および ggplot2というライブラリを読み込んでいます。これらはデータの操作、時系列解析、グラフ作成に役立ちます。
データの前処理:
- flu_seasonは2017年12月16日から2018年1月16日までの日付を生成しています。
- exclude_datesは特定の日付範囲を除外するために使用されています。
データのフィルタリングと計算:
- cdc_state_countsから州が「Maryland」のデータを抽出しています。
- compute_expected関数を使用して期待値を計算しています。keep.components = TRUEは成分を保持するオプションです。
グラフの作成:
- expected_diagnostic関数を使用して、期待値の診断プロットを作成しています。
- このプロットには以下の要素が含まれています:
  - population: 人口の時系列プロット
  - seasonal: 推定された季節効果のプロット
  - trend: 推定されたトレンドの時系列プロット
  - expected: 期待値の時系列プロット
  - residual: 残差の時系列プロット