SIR (Susceptible-Infected-Recovered) モデルで予測する新型コロナウィルス感染の今後の推移

Susceptible-Infected-Recovered (SIR) モデル

Susceptible-Infected-Recovered (SIR) モデルは、 代表的な感染症モデルの１つです。このモデルでは、人口は次の３つのグループに分かれると仮定し、病が伝播しやがて収束するまでのプロセスを記述します。¹

• 感染可能者（Susceptible）：これから感染する可能性のある人
• 感染者（Infected）：現在感染している人
• 除外者（Recovered）：死亡・回復・隔離された人

それぞれの人口における割合を $x, y, z$ とすると、SIRモデルは次の微分方程式で表されます。

\begin{align}
\dot{x} &= - \beta x y \tag{1} \\
\dot{y} &= \beta x y - \gamma y \tag{2} \\
\dot{z} &= \gamma y \tag{3}
\end{align}

ここで、$\dot{x}$ は変数の時間変化（$\mathrm{d}x/\mathrm{d}t$）を表します。
$\beta$ は感染率を意味するパラメータで、これが高いほど人から人へ感染しやすいことを示します。$\gamma$ は、Iが回復・死亡・隔離などの理由によって、単位時間で感染者でなくなる割合です。病気が治りやすい場合 $\gamma$ は高くなる一方、非常に深刻で感染者がすぐに死に至る場合にもやはり $\gamma$ は高くなります。

SIRモデルの基本的な挙動については 別の記事 で紹介しているのでご参照ください。
この記事では、SIRモデルの現実問題への適用に取り組みます。

モデルを実データに当てはめるためには、パラメータ $\beta, \gamma$ を何らかの方法で推定する必要があります。
推定方法の１つとして、微分方程式を解き（$x, y, z$ を $t$の関数として書く）、それと実データにフィットさせる方法があります。Burghes & Borrie (1981) の 7.4 節では、このモデルを解くことはできないと書かれていますが¹、後に Harko et al. (2004) によって明示的な解が提示されています²。 Toda (2020) では、この解と各国の新型コロナウィルスのデータをフィットさせ、非線形最小二乗法によってモデルパラメータを推定しています³。

この記事では、次節で述べるような別の方法でパラメータを推定を試み、これを用いて新型コロナウィルス問題の今後の推移を占います。

用いたコードは前回の記事のものも含めてこちらのGitHubレポジトリに公開しています。

パラメータの推計手法

パラメータの推定式を作るに辺り、少し計算を行います。式(1)と(3)の比をとると、次の微分方程式が得られます。

\frac{\dot{x}}{\dot{z}} = \frac{\mathrm{d}x}{\mathrm{d}z} = -\frac{\beta x}{\gamma}

これを解くと、

\begin{align}
\frac{\mathrm{d}x}{\mathrm{d}z} &= -\frac{\beta x}{\gamma} \\
\frac{\mathrm{d}x}{x} &= -\frac{\beta}{\gamma} \mathrm{d}z \\
\int \frac{\mathrm{d}x}{x} &= - \int \frac{\beta}{\gamma}  \mathrm{d}z \\
\log x &= - \frac{\beta}{\gamma} z + C \\
z &= C' - \rho \log x,\;\;\; \rho \equiv \frac{\gamma}{\beta} \tag{4}
\end{align}

また、(3)式から、
$$
\dot{z} = \frac{\mathrm{d}z}{\mathrm{d}t} = \gamma y \tag{5}
$$
を得られます。

そこで、$x, y, z$ のデータが得られるのであれば、(4) を用いて $\rho=\gamma/\beta$ を推定し、(5) を用いて $\gamma$ を推定することができます。また、$\beta$ は$\gamma/\rho$ により計算することができます。
いずれも線形回帰により推定できる式なので、計算は比較的シンプルに実装できます。

まずは、このアイデアを検証するためシミュレーションデータに当てはめてみます。

library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)

sir_simulate <- function(beta, gamma, maxt=150, y0=1e-6, z0=0) {
  x0 <- 1 - y0 - z0
  x <- numeric(1+maxt)
  y <- numeric(1+maxt)
  z <- numeric(1+maxt)
  x[1] <- x0
  y[1] <- y0
  z[1] <- z0
  for (i in 1:maxt) {
    x[i+1] <- x[i] - beta*x[i]*y[i]
    y[i+1] <- y[i] + beta*x[i]*y[i] - gamma*y[i]
    z[i+1] <- z[i] + gamma*y[i]
  }
  data.frame(t=0:maxt, susceptible=x, infected=y, recovered=z)
}

estimate_sir_params <- function(x, y, z) {
  m1 <- lm(z ~ log(x))
  delta_z <- (lead(z) - lag(z))/2
  m2 <- lm(delta_z ~ y - 1)
  rho <- -m1$coef[2] %>% unname()
  gamma <- m2$coef[1] %>% unname()
  beta <- gamma / rho %>% unname()
  
  # 標準誤差の計算。２つの回帰誤差の無相間を仮定し
  rho_se <- sqrt(diag(vcov(m1)))[2] %>% unname()
  gamma_se <- sqrt(diag(vcov(m2)))[1] %>% unname()
  beta_se <- sqrt(gamma_se^2 / rho^2 + rho_se^2 * (gamma/rho^2)^2) %>% unname()

  out <- list(beta=beta, gamma=gamma, rho=rho, 
              beta_se=beta_se, gamma_se=gamma_se, rho_se=rho_se)
  out
}

d <- sir_simulate(0.3, 0.1)
estimate_sir_params(d$susceptible, d$infected, d$recovered)

結果、悪くない推定値が得られました。$\beta$ の推定値がやや高めになっていますが、おそらく微分方程式で記述されたモデルを差分に置き換えてシミュレーションしている影響ではないかと思います。

$beta
[1] 0.3098709

$gamma
[1] 0.09990712

$rho
[1] 0.3224153

$beta_se
[1] 0.0007799893

$gamma_se
[1] 0.0002503641

$rho_se
[1] 7.638358e-05

国別の感染者推移の予測

データ

推定には、各国のデータを取りまとめてくれている、COVID-19 というGitHubレポジトリを利用します。
このレポジトリでは、国別・日別で「感染確認者」「回復者」「死者」の数が得られます。
モデルとの対応としては、

• 除外者 (R): 「回復者」+「死者」
• 感染者 (I): 「感染確認者」-「回復者」-「死者」
• 感染可能者 (S): それ以外

と仮定します。各国の人口は Wikipedia | List of countries and dependencies by population を参照し、全ての人が S, I, R のいずれかに属すると仮定します。

下記のコードではデータをダウンロードし、人口と結合し推定に便利な形に整形しています。

load_csse_data <- function() {
  y <- read.csv("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv",
                stringsAsFactors=FALSE) %>%
    select(-c(Lat, Long)) %>%
    pivot_longer(-c(Province.State, Country.Region), names_to="date", values_to="confirmed") %>%
    mutate(date=as.Date(date, format="X%m.%d.%y"))
  z1 <- read.csv("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv",
                stringsAsFactors=FALSE) %>%
    select(-c(Lat, Long)) %>%
    pivot_longer(-c(Province.State, Country.Region), names_to="date", values_to="recovered") %>%
    mutate(date=as.Date(date, format="X%m.%d.%y"))
  z2 <- read.csv("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv",
                 stringsAsFactors=FALSE) %>%
    select(-c(Lat, Long)) %>%
    pivot_longer(-c(Province.State, Country.Region), names_to="date", values_to="dead") %>%
    mutate(date=as.Date(date, format="X%m.%d.%y"))

  inner_join(y, z1) %>% inner_join(z2) %>%
    rename(state=Province.State, country=Country.Region)
}
csse <- load_csse_data()

# 一部の国に対象を限定
pop_world <- read.csv(text='country,pop
Italy,60243406
France,67076000
Germany,83149300
Spain,47100396
United Kingdom,66435550
US,329567756
Australia,25667395
Brazil,211345539
Japan,125950000', stringsAsFactors=FALSE)

# 地域別データを国別に集約
tmp <- group_by(csse, date, country) %>%
  summarize(confirmed=sum(confirmed),
            recovered=sum(recovered),
            dead=sum(dead)) 

world_data <- inner_join(pop_world, tmp) %>%
  mutate(infected=confirmed-recovered-dead,
         recovered_or_dead=recovered+dead,
         susceptible=pop-infected-recovered) %>%
  mutate(x=susceptible/pop, y=infected/pop, z=recovered_or_dead/pop)
head(world_data)

下のようなデータが出来上がりました。変数 x, y, z はSIRモデルに出てくる３つのグループの人口比に対応します。

各国の感染者数の推移を見ると、いずれも現在上昇局面にあることが分かります。

ggplot(world_data, aes(date, infected)) +
  geom_line() +
  facet_wrap(vars(country), scales="free_y") +
  theme_bw()

パラメータ推定

国別にパラメータ $\beta, \gamma$ を推定します。各国の政策や事情が日々変化していることから、直近の20日分のデータを用いています。20日が適切な数字かはわかりませんが、結果に大きな差は生じないようです。

world_est <- NULL
for (ctr in unique(world_data$country)) {
  tmp <- filter(world_data, country==ctr) %>% 
    arrange(date) %>%
    tail(20)
  est <- estimate_sir_params(tmp$x, tmp$y, tmp$z)
  tmp <- as.data.frame(est[c("beta", "beta_se", "gamma",
                             "gamma_se", "rho", "rho_se")])
  tmp$country <- ctr
  world_est <- rbind(world_est, tmp)                    
}

ggplot(world_est, aes(beta, gamma)) +
  geom_text(aes(label=country), hjust="inward") +
  theme_bw()

この中では、日本は感染率（$\beta$）・除外率（$\gamma$）ともに低めと推定されました。除外率の高いスペイン・ドイツ・フランスあたりは、死亡者が多い側面もあるかもしれません。

将来予測

推定されたパラメータとSIRモデルを用いて、各国の感染者数（人口比）の推移を予測します。推定誤差の影響を考慮して、下記の5つのシナリオをとっています。

• 「ベースライン」：推定値をそのまま用いる
• 「$\beta$ 高, $\rho$ 高」betaは95％信頼区間の上限値、gammaは上限値を用いる
• 「$\beta$ 高, $\rho$ 低」betaは95％信頼区間の上限値、gammaは下限値を用いる
• 「$\beta$ 低, $\rho$ 高」betaは95％信頼区間の下限値、gammaは上限値を用いる
• 「$\beta$ 低, $\rho$ 低」betaは95％信頼区間の下限値、gammaは下限値を用いる

これらのパラメータを用いて、現時点の最新の $y, z$ を初期値と与えたときの、300日後までの推移をシミュレーションにより算出します。

world_sim <- NULL
for (ctr in unique(world_data$country)) {
  start <- filter(world_data, country==ctr) %>% arrange(date) %>% tail(1)
  est <- filter(world_est, country==ctr)
  tmp1 <- sir_simulate(beta=est$beta, gamma=est$gamma,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="baseline")
  tmp2 <- sir_simulate(beta=est$beta-est$beta_se*2, gamma=est$gamma+est$gamma_se*2,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="low beta/high gamma")
  tmp3 <- sir_simulate(beta=est$beta+est$beta_se*2, gamma=est$gamma-est$gamma_se*2,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="high beta/low gamma")
  tmp4 <- sir_simulate(beta=est$beta+est$beta_se*2, gamma=est$gamma+est$gamma_se*2,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="high beta/high gamma")
  tmp5 <- sir_simulate(beta=est$beta-est$beta_se*2, gamma=est$gamma-est$gamma_se*2,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="low beta/low gamma")
  tmp <- rbind(tmp1, tmp2, tmp3, tmp4, tmp5)

  tmp$country <- ctr
  tmp$date <- start$date + tmp$t
  world_sim <- rbind(world_sim, tmp)                    
}

ggplot(world_data, aes(date, y)) +
  geom_line() +
  geom_line(data=world_sim, aes(date, infected, color=param, linetype=param)) +
  facet_wrap(vars(country), scales="free_y", ncol=3) +
  ylab("Infected") +
  theme_bw() +
  scale_x_date(labels=date_format("%y-%m")) +
  theme(legend.position="bottom", legend.title=element_blank())

イタリア・日本を除く各国はほぼ同様の傾向を見せており、概ね夏頃にピークを迎えて収束に向かうという予測になっています。
イタリアと日本はやや遅めで、秋頃のピークが予測されています。
この主たる要因は $\beta$ の推定値の低さと考えられます。$\beta$（感染率）が低いと感染者が増えるスピードが遅くなるため、集団免疫の獲得するまでに時間がかかるためです。その分除外率（$\gamma$）が高ければ感染者は減っていくのですが、治療薬があるわけでもないため他国と比べて特別低い値にはなっていません。データは検査により確認された感染者数をそのまま用いているため、各国で検査に関わる方針に差がある場合、それが $\beta$ の推定値に影響する可能性があります。

この点は他の国にも言えることで、データとして得られている感染者の数は過小評価されている可能性が高いです。実際にはもっと多くの感染者がいるとするなら、$\beta$ の値もまた過小評価されます。$\beta$ が高くなるとピーク時期と収束までにかかる時間が早まりますから、もう少し早期の収束を期待できるのかもしれません。

日本の感染者推移の予測

データ

同様の分析を日本のデータを用いて行います。全国・東京都・大阪府の3つにわけてパラメータを推定し、今後の動向を予測します。データは、2019-ncov-japan というGitHubレポジトリ利用します。
人口データについては Wikipedia を参照しました。

データとモデルの対応関係としては、下記を仮定します。

• 除外者 (R): 「退院者」+「死者」
• 感染者 (I): 「感染確認者」-「退院者」-「死者」
• 感染可能者 (S): それ以外

jpn <- read.csv("https://raw.githubusercontent.com/swsoyee/2019-ncov-japan/master/Data/detailByRegion.csv",
                stringsAsFactors=FALSE) %>%
  setNames(c("date", "pref", "confirmed", "hospitalized", "recovered", "dead")) %>%
  mutate(date=as.Date(as.character(date), format="%Y%m%d"))

pop_jpn <- read.csv(text="pref,pref_en,pop
全国,Total,127094745
東京都,Tokyo,13515271
大阪府,Osaka,8839469", stringsAsFactors=FALSE)

# 全国合計を計算
tmp <- group_by(jpn, date) %>%
  summarize(confirmed=sum(confirmed),
            hospitalized=sum(hospitalized),
            recovered=sum(recovered),
            dead=sum(dead)) %>%
  ungroup() %>%
  mutate(pref="全国") %>%
  rbind(jpn)

jpn_data <- inner_join(pop_jpn, tmp) %>%
  mutate(infected=confirmed-recovered-dead,
         recovered_or_dead=recovered+dead,
         susceptible=pop-infected-recovered) %>%
  mutate(x=susceptible/pop, y=infected/pop, z=recovered_or_dead/pop) %>%
  mutate(pref_en=factor(pref_en, levels=pop_jpn$pref_en))
head(jpn_data

ggplot(jpn_data, aes(date, infected)) +
  geom_line() +
  facet_wrap(vars(pref_en), scales="free_y") +
  theme_bw()

パラメータ推定

パラメータを推定します。取得可能なデータ期間が短いので全データを用います。

jpn_est <- NULL
for (p in unique(jpn_data$pref_en)) {
  tmp <- filter(jpn_data, pref_en==p) %>% 
    arrange(date) 
  est <- estimate_sir_params(tmp$x, tmp$y, tmp$z)
  tmp <- as.data.frame(est[c("beta", "beta_se", "gamma", "gamma_se", "rho", "rho_se")])
  tmp$pref_en <- p
  jpn_est <- rbind(jpn_est, tmp)  
}

pivot_longer(jpn_est, -pref_en, names_to="variable", values_to="value") %>%
  mutate(type=ifelse(grepl("_se", variable), "se", "estimate")) %>%
  mutate(variable=sub("_se", "", variable)) %>%
  pivot_wider(id_cols=c(pref_en, variable), names_from=type, values_from=value) %>%
  mutate(lb=estimate-2*se, ub=estimate+2*se) %>%
ggplot(aes(pref_en)) +
  geom_segment(aes(y=lb, yend=ub, xend=pref_en)) +
  geom_point(aes(y=estimate)) +
  facet_wrap(vars(variable), scales="free_x") +
  xlab(element_blank()) + ylab("Estimates") +
  coord_flip() +
  theme_bw()

感染率（$\beta$）はほとんど変わりませんが、東京で除外率（$\gamma$, 回復するか死亡する割合）が低いという結果になりました。

将来予測

推定されたパラメータを用いて将来予測を行います。国別の時と同様の5つのシナリオを用いています。

jpn_sim <- NULL
for (p in unique(jpn_data$pref_en)) {
  start <- filter(jpn_data, pref_en==p) %>% arrange(date) %>% tail(1)
  est <- filter(jpn_est, pref_en==p)
  tmp1 <- sir_simulate(beta=est$beta, gamma=est$gamma,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="baseline")
  tmp2 <- sir_simulate(beta=est$beta-est$beta_se*2, gamma=est$gamma+est$gamma_se*2,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="low beta/high gamma")
  tmp3 <- sir_simulate(beta=est$beta+est$beta_se*2, gamma=est$gamma-est$gamma_se*2,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="high beta/low gamma")
  tmp4 <- sir_simulate(beta=est$beta+est$beta_se*2, gamma=est$gamma+est$gamma_se*2,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="high beta/high gamma")
  tmp5 <- sir_simulate(beta=est$beta-est$beta_se*2, gamma=est$gamma-est$gamma_se*2,
                       y0=start$y, z0=start$z, maxt=300) %>%
    mutate(param="low beta/low gamma")
  tmp <- rbind(tmp1, tmp2, tmp3, tmp4, tmp5)

  tmp$pref_en <- p
  tmp$date <- start$date + tmp$t
  jpn_sim <- rbind(jpn_sim, tmp)                    
}
jpn_sim$pref_en <- factor(jpn_sim$pref_en, levels=levels(jpn_data$pref_en))

options(repr.plot.width=7.5, repr.plot.height=3)
ggplot(jpn_data, aes(date, y)) +
  geom_line() +
  geom_line(data=jpn_sim, aes(date, infected, color=param, linetype=param)) +
  facet_wrap(vars(pref_en), scales="free_y", ncol=3) +
  ylab("Infected") +
  theme_bw() +
  scale_x_date(labels=date_format("%y-%m")) +
  theme(legend.position="bottom", legend.title=element_blank())

国別の結果と似通った結果になっており、ベースラインでは概ね夏の終わりから秋にかけてピークを迎えると予測されています。$\beta$ を高めに設定したシナリオでは夏頃のピーク、$\beta$ を低めに想定したシナリオでは、ゆるやかな感染が長期化するという結果になっています。
もしデータが感染者の数を捕捉しきれておらず、$\beta$ の推定値が過小評価されているとするならば、予測結果よりもう少し早い収束も期待できるかもしれません。

結び

当然のことながら、この分析はあくまで１つの試みで、現実はまったく違う動きをする可能性が十分ありますし、実際そうなる確率の方が高いでしょう。
SIRモデルを仮定していますが、これが新型コロナウウィルスの伝染に適したモデルであるかはわかりません。また、各国の検査・症例報告にもとづくデータを用いており、これは現時点で手に入るベストに近いものだとは思うものの、必ずしもモデルの想定する変数に対応していないところもあります。

仮にこの分析手法が正しいとしても考慮に入れていない点もあります。１つは各国の政策変化です。予測にあたっては現時点のパラメータ推定値をあてはめていますが、実際には規制等の実施によりパラメータは日々変化していきます。また、ワクチン・薬品・治療法などの医療技術の開発・進歩も想定されていません。SIRモデルでの収束には、社会の大半が感染を経験することによる集団免疫の獲得が必要となります。ですが、理想的にはそれよりも先に医学的な進歩によって感染率が劇的に下がったり治癒率があがったりするかもしれません。これもモデルが組み入れていないパラメータの変化の一種です。

1. Burghes & Borrie (1981) Modelling with Differential Equations（垣田・大町訳 (1990)『微分方程式で数学モデルを作ろう』）の7.4節を参照 ↩ ↩²

2. Harko T., Lobo, S.N.F, Mak, M.K. (2004) Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates, Applied Mathematics and Computation 236. ↩

3. Toda, A.A. (2020) Susceptible-Infected-Recovered (SIR) dynamics of COVID-19 and Economic Impact, arXiv:2003.11221. ↩