24 エピデミックモデリング
24.1 概要
流行のモデリングのためのツールが増加してきており、最小の労力でかなり複雑な分析を行うことができるようになっています。このセクションでは、これらのツールをどのように使用するかについての概要を説明します：
- 実効再生産数（effective reproduction number）Rt や倍化時間（doubling time）などに関連する統計量を推定する
- 将来のインシデンスの短期的な予測（プロジェクション）を実施する
ここではこれらの基礎となっている方法論や統計学的な手法の概要を説明する意図はありませんので、関連論文へのリンクについては Resources tab を参照してください。ここで説明するツールを使用する前に、手法を理解しておくことで、正確に結果を解釈できるようになります。
以下はこのセクションで作成するアウトプットの１つの例です。
24.2 準備
Rt の推定には EpiNow パッケージと EpiEstim パッケージという2つの異なる手法を使用し、症例（ケース）の発生数の予測には projections パッケージを用います。以下のコードチャンクは、（本章の）分析に必要なパッケージのローディングを示しています。
以下のコードを実行すると、分析に必要なパッケージが読み込まれます。このハンドブックでは、パッケージを読み込むために、pacman パッケージの p_load() を主に使用しています。p_load() は、必要に応じてパッケージをインストールし、現在の R セッションで使用するためにパッケージを読み込む関数です。また、すでにインストールされたパッケージは、R の基本パッケージである base （以下、base R）の library() を使用して読み込むこともできます。R のパッケージに関する詳細は R の基礎 の章をご覧ください。
pacman::p_load(
rio, # ファイルをインポート
here, # ファイルロケーター
tidyverse, # データマネジメント + ggplot2 のグラフィックス
epicontacts, # トランスミッションネットワークの分析
EpiNow2, # Rt 推定
EpiEstim, # Rt 推定
projections, # 発生数のプロジェクション
incidence2, # 発生データの取り扱い
epitrix, # 便利な epi の機能
distcrete # 離散的な遅れの分布
)エボラ出血熱の流行をシミュレートしたデータセットをインポートします。お手元の環境でこの章の内容を実行したい方は、 クリックして「前処理された」ラインリスト（linelist）データをダウンロードしてください>（.rds 形式で取得できます）。データは rio パッケージの import() を利用してインポートしましょう（rio パッケージは、.xlsx、.csv、.rds など様々な種類のファイルを取り扱うことができます。詳細は、インポートとエクスポート の章をご覧ください）。
# クリーンなラインリストの取り込み
linelist <- import("linelist_cleaned.rds")24.3 Rt の推定
EpiNow2 vs. EpiEstim
再生産数 R は、疾病の感染性を示す指標であり、感染者 1 人あたりの二次感染者数の期待値として定義されます。感受性保持者しかいないような集団では、この値は基本再生産数 R0 を意味します。しかし、集団内の感受性保持者はアウトブレイクやパンデミックの期間中に変化し、また様々な対策が講じられるため、感染性の指標として最も一般的に使用されるのは実効再生産数です（ある時刻 t における感染者1人あたりの二次感染者数の期待値）。
EpiNow2 パッケージは最も洗練された Rt 推定のためのフレームワークを提供しています。ほかに一般的に使用されている EpiEstim パッケージと比較して、2つの重要な利点があります：
報告の遅れを考慮しているので、直近のデータが不完全の場合であっても Rt を推定することができます。
報告日ではなく、感染日に基づいて Rt を推定するので、遅れを生じずすぐに Rt の変化として介入の効果が反映されま。
しかし、2つの重要なデメリットもあります：
- 世代時間の分布（一次感染者から二次感染者までの遅れの分布）、潜伏期間の分布（感染から症状発現までの遅れの分布）、およびデータに関連するその他の遅れの分布（例えば、報告の日付がある場合、症状発現から報告までの遅れの分布）に関する知識が必要です。これにより、より正確な Rt を推定できますが、EpiEstim パッケージは発症間隔発症間隔（一次感染者の症状発現から二次感染者の症状発現までの遅れの分布）のみを必要とし、これのみが唯一利用可能な分布である場合があります。
- EpiNow2 パッケージは EpiEstim パッケージに比べて著しく遅く、約100～1000倍の差があるといわれています！例えば、このセクションで利用するサンプルアウトブレイクにおける Rt 推定には、約4時間かかります（これは高い精度を確保するために多数の反復処理を行ったためであり、必要に応じて短縮することも可能ですが、アルゴリズムが一般的に遅いという点は変わりません）。定期的に Rt の推定値を更新している場合は、この方法は現実的ではないかもしれません。
そのため、どのパッケージを選ぶかは、利用できるデータや時間、計算資源によります。
EpiNow2
遅れの分布の推定
EpiNow2 パッケージの実行に必要な遅れの分布は、手持ちのデータによって異なります。基本的には、感染日から Rt 推定に使用したいイベント日までの遅れを記述できるものである必要があります。もし発症日を使っている場合、これは単に潜伏期間の分布となります。報告日を使用している場合は、感染から報告までの遅れの分布が必要です。この分布はなかなか直接知ることができないため、EpiNow2 パッケージでは複数の遅れの分布をつなぎ合わせることができます。この場合、感染から症状発現までの遅れ（例えば潜伏期間、これは既知であることが多いです）と、症状発現から報告までの遅れ（これは自分でデータから推定できる場合が多いです）です。
例のラインリストではすべての症例について発症日がわかっているので、データ（例えば症状の発現日など）を感染日に結びつけるためには、潜伏期間の分布が必要になります。この分布は、データから推定するか、既存文献から値を引用することができます。
エボラ出血熱の潜伏期間を平均9.1日、標準偏差7.3日、最大値を30日とする文献からの推定値（引用論文）は、以下のように規定されます：
incubation_period_lit <- list(
mean = log(9.1),
mean_sd = log(0.1),
sd = log(7.3),
sd_sd = log(0.1),
max = 30
)EpiNow2 パッケージでは、これらの遅れの分布が対数（log）スケールで提供される必要があり、そのため、各値に log がついていることに注意してください（紛らわしいことに、自然スケールで提供されなければならない max パラメータを除く）。mean_sd と sd_sd は、平均値と標準偏差の推定値の標準偏差を定義します。上記のケースではこれらは知られていないため、かなり恣意的な値である0.1を選択しました。
今回のｂ分析ではその代わりに、bootstrapped_dist_fit() を用いて、ラインリストから潜伏期間の分布を推定しました。
## 潜伏期間の推定
incubation_period <- bootstrapped_dist_fit(
linelist$date_onset - linelist$date_infection,
dist = "lognormal",
max_value = 100,
bootstraps = 1
)もう一つ必要な分布は、世代時間です。感染時刻と感染伝播のリンクに関するデータがあるので、感染者と被感染者のペアの感染時刻の遅れを計算することで、ラインリストからこの分布を推定することができます。これには epicontacts パッケージにある便利な get_pairwise() を使います。この関数を使うと、感染ペア間のラインリスト上の2組の特性の違いを計算することができます。epicontacts オブジェクトを作成します（詳しくは 感染連鎖の図式化 の章を参照してください）：
## コンタクトの作成
contacts <- linelist %>%
transmute(
from = infector,
to = case_id
) %>%
drop_na()
## epicontacts オブジェクトの作成
epic <- make_epicontacts(
linelist = linelist,
contacts = contacts,
directed = TRUE
)次に、get_pairwise で計算した感染ペア間の感染時刻の差をガンマ分布にあてはめました：
## ガンマ分布に従う世代時間の推定
generation_time <- bootstrapped_dist_fit(
get_pairwise(epic, "date_infection"),
dist = "gamma",
max_value = 20,
bootstraps = 1
)EpiNow2 パッケージの実行
あとはラインリストから日々のインシデンスを計算するだけですが、dplyr パッケージの group_by() と n() で簡単にできます。EpiNow2 パッケージでは、列名が date と confirm でなければならないことに注意してください。
そして、epinow() を使って Rt を推定することができます。入力に関して、いくつかの注意点を挙げます：
-
delaysの引数には、任意の数の「連鎖した」遅れの分布を与えることができます。delay_opts()内でincubation_periodオブジェクトと一緒に入れるだけです。 -
return_outputは、出力ファイルに保存されるのではなく、R の中で貸せされるようになっています。 -
verboseは、進捗状況の読み上げを指定します。 -
horizonは、将来のインシデンスを何日分予測するかを示します。 -
stanの因数に追加のオプションを渡して、推定を実行する期間を指定します。samples数とchains数を増やすと、不確実性の特徴をよく表したより正確な推定値が得られますが、実行には時間がかかります。
## epinow を走らせる
epinow_res <- epinow(
reported_cases = cases,
generation_time = generation_time,
delays = delay_opts(incubation_period),
return_output = TRUE,
verbose = TRUE,
horizon = 21,
stan = stan_opts(samples = 750, chains = 4)
)アウトプットの分析
コードの実行が終了すると、以下のように簡単にサマリーをプロットすることができます。画像をスクロールすると、全体を見ることができます。
## サマリーフィギュアのプロット
plot(epinow_res)
また、様々なサマリー統計量を見ることもできます：
## サマリーテーブル
epinow_res$summary## measure estimate
## 1: New confirmed cases by infection date 4 (2 -- 6)
## 2: Expected change in daily cases Unsure
## 3: Effective reproduction no. 0.88 (0.73 -- 1.1)
## 4: Rate of growth -0.012 (-0.028 -- 0.0052)
## 5: Doubling/halving time (days) -60 (130 -- -25)
## numeric_estimate
## 1: <data.table[1x9]>
## 2: 0.56
## 3: <data.table[1x9]>
## 4: <data.table[1x9]>
## 5: <data.table[1x9]>さらなる分析やカスタムプロットのために $estimates$summarised を介して要約された毎日の推定値にアクセスすることができます。これをデフォルトの data.table から、dplyr パッケージで使いやすいように tibble に変換します。
## サマリーを抽出して、tibble に変換
estimates <- as_tibble(epinow_res$estimates$summarised)
estimates例として、倍化時間と Rt をプロットしてみましょう。極端に高い倍化時間をプロットしないように、Rt が1を大きく上回っている流行の最初の数か月だけを見ています。
log(2)/growth_rate という計算式を用いて、推定された成長率（growth rate）から倍化時間を算出しています。
## 中央値プロットのために横型のデータフレームを作ります
df_wide <- estimates %>%
filter(
variable %in% c("growth_rate", "R"),
date < as.Date("2014-09-01")
) %>%
## 成長率を倍化時間に変換
mutate(
across(
c(median, lower_90:upper_90),
~ case_when(
variable == "growth_rate" ~ log(2)/.x,
TRUE ~ .x
)
),
## 変形を反映した変数名の変更
variable = replace(variable, variable == "growth_rate", "doubling_time")
)
## 分位値プロットのために縦長のデータフレームを作る
df_long <- df_wide %>%
## ここでは、マッチした分位値を利用します（例：lower_90 から upper_90）
pivot_longer(
lower_90:upper_90,
names_to = c(".value", "quantile"),
names_pattern = "(.+)_(.+)"
)
## プロットする
ggplot() +
geom_ribbon(
data = df_long,
aes(x = date, ymin = lower, ymax = upper, alpha = quantile),
color = NA
) +
geom_line(
data = df_wide,
aes(x = date, y = median)
) +
## label_parsedを使用して、添え字ラベルを許可する
facet_wrap(
~ variable,
ncol = 1,
scales = "free_y",
labeller = as_labeller(c(R = "R[t]", doubling_time = "Doubling~time"), label_parsed),
strip.position = 'left'
) +
## 分位値の透明度を手動で定義する
scale_alpha_manual(
values = c(`20` = 0.7, `50` = 0.4, `90` = 0.2),
labels = function(x) paste0(x, "%")
) +
labs(
x = NULL,
y = NULL,
alpha = "Credibel\ninterval"
) +
scale_x_date(
date_breaks = "1 month",
date_labels = "%b %d\n%Y"
) +
theme_minimal(base_size = 14) +
theme(
strip.background = element_blank(),
strip.placement = 'outside'
)
EpiEstim
EpiEstim パッケージを走らせるために、日々のインシデンスのデータを提供し、発症間隔（一次症例と二次症例の症状発現までの遅れの分布）を指定する必要があります。
インシデンスデータは、ベクトル、データフレーム、またはオリジナルの incidence パッケージから得られた incidence オブジェクトとして、EpiEstim パッケージに提供できます。輸入感染例とローカルで感染した例を区別することもできます：詳細は ?estimate_R のドキュメントを参照してください。
ここでは incidence2 パッケージを使って、インプットを作成します。incidence2 パッケージの例については 流行曲線（エピカーブ） の章を参照してください。incidence2 パッケージには estimateR() が期待するインプットとは完全に一致しないアップデートがあったため、いくつかの小さな追加手順が必要となります。incidence オブジェクトは日付とそれぞれのケースカウントをもつ tibble で構成されています。tidyr パッケージの complete() を使用して、すべての日付が含まれていることを確認し（症例がない日も含む）、後のステップで estimate_R() で期待されるものと一致するように列を rename() します。
## 発症日からインシデンスを得る
cases <- incidence2::incidence(linelist, date_index = "date_onset") %>% # 日ごとにケースカウントを得る
tidyr::complete(date_index = seq.Date( # すべての日付が表示されていることを確認
from = min(date_index, na.rm = T),
to = max(date_index, na.rm=T),
by = "day"),
fill = list(count = 0)) %>% # NA カウントを0に変換する
rename(I = count, # estimateRで期待されるな目に変更
dates = date_index)このパッケージにはは発症間隔を指定するためのいくつかのオプションがあり、その詳細はドキュメントの ?estimate_R に記載されています。ここではそのうちの2つを取り上げます。
文献から引用した発症間隔の推定値の使用
オプションの method = "parametric_si" を使用すると、 make_config() で作成した config オブジェクトに発症間隔の平均値と標準偏差を手動で指定することができます。ここでは、この論文で定義されている平均値12.0、標準偏差5.2を使用しています。
## config の作成
config_lit <- make_config(
mean_si = 12.0,
std_si = 5.2
)そして、estimate_R() で Rt の推定をすることができます：
cases <- cases %>%
filter(!is.na(date))
#create a dataframe for the function estimate_R()
cases_incidence <- data.frame(dates = seq.Date(from = min(cases$dates),
to = max(cases$dates),
by = 1))
cases_incidence <- left_join(cases_incidence, cases) %>%
select(dates, I) %>%
mutate(I = ifelse(is.na(I), 0, I))## Joining with `by = join_by(dates)`
epiestim_res_lit <- estimate_R(
incid = cases_incidence,
method = "parametric_si",
config = config_lit
)## Default config will estimate R on weekly sliding windows.
## To change this change the t_start and t_end arguments.あとはアウトプットのサマリーをプロットします：
plot(epiestim_res_lit)
データから推定した発症間隔の推定値の使用
症状の発症日と感染伝播のリンクデータがあるので、感染者と被感染者のペアの発症日の遅れを計算することで、ラインリストから発症間隔を推定することもできます。EpiNow2 パッケージのセクションで行ったように、epicontacts パッケージの get_pairwise() を使います。この関数は感染ペア間のラインリスト上の2組の特性の違いを計算することができます。まず epicontacts オブジェクトを作成します（詳細は 感染連鎖の図式化 の章を参照）：
## コンタクトの作成
contacts <- linelist %>%
transmute(
from = infector,
to = case_id
) %>%
drop_na()
## epicontacts オブジェクトの作成
epic <- make_epicontacts(
linelist = linelist,
contacts = contacts,
directed = TRUE
)次に、get_pairwise() を用いて感染ペア間の発症日の差をガンマ分布に当てはめます。離散化された分布が必要とされるため、このフィッティング手順には epitrix パッケージの便利な fit_disc_gamma() を使用します。
## ガンマ分布に従う発症間隔の推定
serial_interval <- fit_disc_gamma(get_pairwise(epic, "date_onset"))その後に、config オブジェクトに情報を与え、 EpiEstim を再実行して結果を描画しましょう。
## config の作成
config_emp <- make_config(
mean_si = serial_interval$mu,
std_si = serial_interval$sd
)
## epiestim を走らせる
epiestim_res_emp <- estimate_R(
incid = cases_incidence,
method = "parametric_si",
config = config_emp
)## Default config will estimate R on weekly sliding windows.
## To change this change the t_start and t_end arguments.
## アウトプットをプロットする
plot(epiestim_res_emp)
推定時間枠（ウィンドウ）の設定
デフォルトのオプションでは週単位の平滑化された推定値を提供していますが、正確な推定値を得るために、アウトブレイク初期に Rt を推定していることを警告する機能もあります。以下に示すように、推定の開始日を遅く設定することで、変更できます。残念ながら、EpiEstim パッケージはこれらの推定時間を指定するのに非常に面倒な方法しか提供しておらず、そのためには各時間ウィンドウの開始日と終了日を参照する整数のベクトルを提供しなければなりません。
## 6月1日から始まる日付のベクトルを定義する
start_dates <- seq.Date(
as.Date("2014-06-01"),
max(cases$dates) - 7,
by = 1
) %>%
## 数字型に変換するために開始日を引く
`-`(min(cases$dates)) %>%
## convert to integer
as.integer()
## 1週間の平滑化ウィンドウに6日分を追加する
end_dates <- start_dates + 6
## config を作成する
config_partial <- make_config(
mean_si = 12.0,
std_si = 5.2,
t_start = start_dates,
t_end = end_dates
)ここで、EpiEstim を再び実行してみると、推定値は6月からしか出ないことがわかります：
## epiestim を走らせる
epiestim_res_partial <- estimate_R(
incid = cases_incidence,
method = "parametric_si",
config = config_partial
)
## アウトプットをプロットする
plot(epiestim_res_partial)
アウトプットの分析
主なアウトプットは $R でアクセスできます。Rt のプロットと Rt とその日に報告された症例数で与えられた「伝播能力」の指標を作成します（これは次世代の感染者数の期待値として表されます）。
## 中央値のために横型のデータフレームを作成します
df_wide <- epiestim_res_lit$R %>%
rename_all(clean_labels) %>%
rename(
lower_95_r = quantile_0_025_r,
lower_90_r = quantile_0_05_r,
lower_50_r = quantile_0_25_r,
upper_50_r = quantile_0_75_r,
upper_90_r = quantile_0_95_r,
upper_95_r = quantile_0_975_r,
) %>%
mutate(
## t_startからt_endまでの日付の中央値を抽出する
dates = epiestim_res_emp$dates[round(map2_dbl(t_start, t_end, median))],
var = "R[t]"
) %>%
## 日々の発生データを統合する
left_join(cases, "dates") %>%
## すべてのr推定値のリスクを計算する
mutate(
across(
lower_95_r:upper_95_r,
~ .x*I,
.names = "{str_replace(.col, '_r', '_risk')}"
)
) %>%
## r推定値とリスク推定値を分離する
pivot_longer(
contains("median"),
names_to = c(".value", "variable"),
names_pattern = "(.+)_(.+)"
) %>%
## 因子（ファクタ）のレベルを割り当てる
mutate(variable = factor(variable, c("risk", "r")))
## クォンタイル（分位値）から縦型のデータフレームを作成する
df_long <- df_wide %>%
select(-variable, -median) %>%
## r/riskの推定値と分位値を分離する
pivot_longer(
contains(c("lower", "upper")),
names_to = c(".value", "quantile", "variable"),
names_pattern = "(.+)_(.+)_(.+)"
) %>%
mutate(variable = factor(variable, c("risk", "r")))
## プロットを作成する
ggplot() +
geom_ribbon(
data = df_long,
aes(x = dates, ymin = lower, ymax = upper, alpha = quantile),
color = NA
) +
geom_line(
data = df_wide,
aes(x = dates, y = median),
alpha = 0.2
) +
## label_parsed を使用して、添え字ラベルをつける
facet_wrap(
~ variable,
ncol = 1,
scales = "free_y",
labeller = as_labeller(c(r = "R[t]", risk = "Transmission~potential"), label_parsed),
strip.position = 'left'
) +
## 分位値の透明度を手動で定義する
scale_alpha_manual(
values = c(`50` = 0.7, `90` = 0.4, `95` = 0.2),
labels = function(x) paste0(x, "%")
) +
labs(
x = NULL,
y = NULL,
alpha = "Credible\ninterval"
) +
scale_x_date(
date_breaks = "1 month",
date_labels = "%b %d\n%Y"
) +
theme_minimal(base_size = 14) +
theme(
strip.background = element_blank(),
strip.placement = 'outside'
)
24.4 インシデンス（発生数）の予測（プロジェクション）
EpiNow2
Rt の推定に加えて、EpiNow2 パッケージは EpiSoon パッケージとの統合により、Rt のプロジェクションや症例数のプロジェクションもサポートします。必要なのは、epinow() の呼び出しで horizon 引数を指定して、何日先までプロジェクションしたいかを示すことだけです。このセクションでは、epinow_res オブジェクトに格納されている分析のアウトプットをプロットします。
## プロットの一番小さい日付を設定する
min_date <- as.Date("2015-03-01")
## 要約された推定値を抽出する
estimates <- as_tibble(epinow_res$estimates$summarised)
## 発生症例数の生データを抽出する
observations <- as_tibble(epinow_res$estimates$observations) %>%
filter(date > min_date)
## 症例数の予測値を抽出する
df_wide <- estimates %>%
filter(
variable == "reported_cases",
type == "forecast",
date > min_date
)
## 分位値プロットのためにさらに横型のフォーマットに変換する
df_long <- df_wide %>%
## ここで、分位値を一致させる（たとえば lower_90 から upper_90）
pivot_longer(
lower_90:upper_90,
names_to = c(".value", "quantile"),
names_pattern = "(.+)_(.+)"
)
## プロットする
ggplot() +
geom_histogram(
data = observations,
aes(x = date, y = confirm),
stat = 'identity',
binwidth = 1
) +
geom_ribbon(
data = df_long,
aes(x = date, ymin = lower, ymax = upper, alpha = quantile),
color = NA
) +
geom_line(
data = df_wide,
aes(x = date, y = median)
) +
geom_vline(xintercept = min(df_long$date), linetype = 2) +
## 分位値の透明度を手動で定義する
scale_alpha_manual(
values = c(`20` = 0.7, `50` = 0.4, `90` = 0.2),
labels = function(x) paste0(x, "%")
) +
labs(
x = NULL,
y = "Daily reported cases",
alpha = "Credible\ninterval"
) +
scale_x_date(
date_breaks = "1 month",
date_labels = "%b %d\n%Y"
) +
theme_minimal(base_size = 14)
予測（プロジェクション）
RECON が開発した projections パッケージでは、実効再生産数 Rt と発症間隔の知識だけで非常に簡単に短期的なインシデンスの予測を行うことができます。ここでは、文献から得られた発症間隔の推定値を使用する方法と、ラインリストから得られた独自の推定値を使用する方法について説明します。
文献での発症間隔の推定値を利用
projections パッケージには、distcrete パッケージに含まれる distcrete クラスの離散化された発症間隔の分布が必要です。ここでは、この論文 で定義された平均12.0、標準偏差5.2のガンマ分布を使用します。これらの値をガンマ分布に必要な shape（形状）および scale（尺度）パラメータに変換するために、epitrix パッケージの gamma_mucv2shapescale() を使用します。
## 変動係数から形状と尺度パラメータを得ることができます
##（例：平均値に対する標準偏差の比）
shapescale <- epitrix::gamma_mucv2shapescale(mu = 12.0, cv = 5.2/12)
## distcrete オブジェクトを作成する
serial_interval_lit <- distcrete::distcrete(
name = "gamma",
interval = 1,
shape = shapescale$shape,
scale = shapescale$scale
)ここでは、発症間隔が正しいことを確認するために簡単なチェックを行います。先ほど定義したガンマ分布の確率密度に $d でアクセスしますが、これは dgamma を呼び出した時と同じです。
## 発症間隔が正しいことを確認する
qplot(
x = 0:50, y = serial_interval_lit$d(0:50), geom = "area",
xlab = "Serial interval", ylab = "Density"
)## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
データによる発症間隔の推定値を利用
我々は症状の発症日と感染リンク（transmission links）のデータがあるので、感染者と被感染者のペアの発症日の遅れを計算することで、ラインリストから発症間隔を推定することもできます。EpiNow2 のセクションで行ったように、epicontacts パッケージの get_pairwise() を用います。この関数によって、感染ペア（transmission pairs）間のラインリスト上の2組の特性の違いを計算することができます。まず epicontacts オブジェクトを作成します（感染連鎖の図式化 の章を参照）：
## コンタクトの作成
contacts <- linelist %>%
transmute(
from = infector,
to = case_id
) %>%
drop_na()
## epicontacts オブジェクトの作成
epic <- make_epicontacts(
linelist = linelist,
contacts = contacts,
directed = TRUE
)次に、get_pairwise() を用いて、計算した感染ペアの発症日の差をガンマ分布にあてはめます。離散化された分布が必要なため、この適合（fitting）手順には epitrix パッケージの便利な fit_disc_gamma() を使います。
## ガンマ分布に従う発症間隔の分布を推定する
serial_interval <- fit_disc_gamma(get_pairwise(epic, "date_onset"))
## 推定値
serial_interval[c("mu", "sd")]## $mu
## [1] 11.51047
##
## $sd
## [1] 7.696056発生数（インシデンス）の予測（プロジェクション）
将来のインシデンスをプロジェクションするためには、過去の発生数を incidence オブジェクトの形で値供することに加えて、妥当な Rt 値のサンプルを与える必要があります。前のセクション（Rt推定）で EpiEstim によって作成され、epiestim_res_emp オブジェクトに格納された Rt の推定値を使用して、インシデンスの予測値を作成します。以下のコードでは、アウトブレイクの最後の時間ウィンドウの Rt の平均値と標準偏差の推定値を抜き出し（ベクトルの最後の要素にアクセスするために tail() を使用）、rgamma()を使用してガンマ関数から10000の値をシミュレーションします。また、事前予測に使用したいRt 値の独自のベクトルを提供することもできます。
## 発症日からincidenceオブジェクトを作成
inc <- incidence::incidence(linelist$date_onset)## 256 missing observations were removed.
## 最新の推定値から妥当な r 値を抽出する
mean_r <- tail(epiestim_res_emp$R$`Mean(R)`, 1)
sd_r <- tail(epiestim_res_emp$R$`Std(R)`, 1)
shapescale <- gamma_mucv2shapescale(mu = mean_r, cv = sd_r/mean_r)
plausible_r <- rgamma(1000, shape = shapescale$shape, scale = shapescale$scale)
## 分布をチェック
qplot(x = plausible_r, geom = "histogram", xlab = expression(R[t]), ylab = "Counts")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
そして、project() を使って、実際の予測を行います。n_days 引数で何日間プロジェクションするかを指定し、n_sim 引数でシミュレーション回数を指定します。
## プロジェクションの作成
proj <- project(
x = inc,
R = plausible_r,
si = serial_interval$distribution,
n_days = 21,
n_sim = 1000
)そして、plot() と add_projections() を使って、インシデンスとプロジェクションを簡単にプロットすることができます。各括弧演算子（[]）を使用すると、incidenceオブジェクトを簡単に分けることができ、最近のケースのみを表示することができます。
また、アウトプットをデータフレームに変換することで、日々の症例数の生（raw）の推定値を簡単に取り出すことができます。
## 生データをデータフレームに変換する
proj_df <- as.data.frame(proj)
proj_df