ＲでCOVID-19の感染者について予測する

もうｎ番煎じという気もしますが今回はRを用いて今世間をにぎわせているCOVID-19の感染者予測やシミュレーションを行ってみたいと思います。

データの準備

まず最初にデータを準備したいと思います。
今回は用いるデータはこちらからダウンロードできるデータを用いました。
次にデータをRに読み込みたいと思います。以下では2020年4月21日現在のデータを用いて解析しました。このデータは日々更新されています。

library(readr)
cases <- read_csv("COVID-19.csv")

データの可視化

次にデータ分析の最も基本ともいえるデータの可視化を行いたいと思います。とはいえ、こんなにもデータのカテゴリーがあれば、何から図示していけばよいのかわからないので、とりあえず感染者数が多いとされている東京都、神奈川県、大阪府について感染者数をグラフにしてみます。

x<-c(length(which(cases$受診都道府県=="東京都")),length(which(cases$受診都道府県=="大阪府")),length(which(cases$受診都道府県=="神奈川県")))
names(x)<-c("東京都","大阪府","神奈川県")
barplot(x,xlab="prefecture",ylab="population")

これにより以下のようなグラフが得られました。

次に年代別の感染者も図示してみましょう。(若者の感染者が目立つといわれていたので)

a<-length(which(cases$年代<=20))
b<-length(which((cases$年代>20)&(cases$年代<=30)))
c<-length(which((cases$年代>30)&(cases$年代<=40)))
d<-length(which((cases$年代>40)&(cases$年代<=50)))
e<-length(which((cases$年代>50)&(cases$年代<=60)))
f<-length(which((cases$年代>60)&(cases$年代<=70)))
g<-length(which((cases$年代>70)))
data<-c(a,b,c,d,e,f,g)
barplot=barplot(data,names.arg=c("<=20","<=30","<=40","<=50",'<=60','<=70','70<'),ylim=c(0,2500))

うーむ…
もちろん年代別の人口も異なるのでこのグラフだけからでは何とも言えない…
そこで次にデータの可視化とは趣旨が異なりますが、実際の日本の人口比率と感染者の年齢比率の差の検定を行ってみたいと思います。
実際の日本の人口比率はこちらより拝借いたしました。(総務所統計局のデータ見てるだけで時間が溶けますね…楽しい…)
Ｒを用いて比率の差の検定を行うやり方はこのページを参考にしました。

prop.test(c(2153,2265),c(12656,11125),correct=T) #20歳以下の比率の差の検定
prop.test(c(1249,1709),c(12656,11125),correct=T) #30歳以下20歳より上での比率の差の検定
prop.test(c(1487,1873),c(12656,11125),correct=T) #40歳以下30歳より上での比率の差の検定
prop.test(c(1582,1287),c(12656,11125),correct=T) #60歳以下50歳より上での比率の差の検定
prop.test(c(1743,1002),c(12656,11125),correct=T) #70歳以下60歳より上での比率の差の検定
prop.test(c(2558,1064),c(12656,11125),correct=T) #70歳以上の比率の差の検定

このようにして得られた結果が以下です.

#20歳以下の比率の差の検定
##  ##  2-sample test for equality of proportions with continuity correction 
##  ## data:  c(2153, 2265) out of c(12656, 11125) 
## X-squared = 43.648, df = 1, p-value = 3.93e-11 
## alternative hypothesis: two.sided ## 95 percent confidence interval: 
##  -0.04350483 -0.02345230 
## sample estimates: 
##    prop 1    prop 2  
## 0.1701169 0.2035955 

#30歳以下 20歳より上での比率の差の検定
 
##  ##  2-sample test for equality of proportions with continuity correction 
##  ## data:  c(1249, 1709) out of c(12656, 11125) 
## X-squared = 163.52, df = 1, p-value < 2.2e-16 
## alternative hypothesis: two.sided ## 95 percent confidence interval: 
##  -0.06349311 -0.04636611 
## sample estimates: 
##     prop 1     prop 2  
## 0.09868837 0.15361798 

#40歳以下 30歳より上での比率の差の検定
 
##  ##  2-sample test for equality of proportions with continuity correction 
##  ## data:  c(1487, 1873) out of c(12656, 11125) 
## X-squared = 125.84, df = 1, p-value < 2.2e-16 
## alternative hypothesis: two.sided ## 95 percent confidence interval: 
##  -0.05988450 -0.04184724 
## sample estimates: 
##    prop 1    prop 2  
## 0.1174937 0.1683596
 

それぞれの群においてp値が0.05を切っていて、たしかに若者(40歳以下)が日本の人口比率に比べて有意に多く感染していることが分かりました。(もちろん、これらからだけでは判断できないですが。近年(とはいってもだいぶ昔からですが)、p値が再現性の低い値であることが話題になっていますね。)

次に50歳以上についても同様なことをやってみましょう。

#60歳以下 50歳より上での比率の差の検定
 
##  ##  2-sample test for equality of proportions with continuity correction 
##  ## data:  c(1582, 1287) out of c(12656, 11125) 
## X-squared = 4.7547, df = 1, p-value = 0.02922 
## alternative hypothesis: two.sided ## 95 percent confidence interval: 
##  0.0009522702 0.0176769433 
## sample estimates: 
##    prop 1    prop 2  
## 0.1250000 0.1156854 

#70歳以下 60歳より上での比率の差の検定
 
##  ##  2-sample test for equality of proportions with continuity correction 
##  ## data:  c(1743, 1002) out of c(12656, 11125) 
## X-squared = 131.21, df = 1, p-value < 2.2e-16 
## alternative hypothesis: two.sided ## 95 percent confidence interval: 
##  0.03954788 0.05575977 
## sample estimates: 
##     prop 1     prop 2  
## 0.13772124 0.09006742 

#70歳以上の比率の差の検定
 
##  ##  2-sample test for equality of proportions with continuity correction 
##  ## data:  c(2558, 1064) out of c(12656, 11125) 
## X-squared = 519.08, df = 1, p-value < 2.2e-16 
## alternative hypothesis: two.sided ## 95 percent confidence interval: 
##  0.09751489 0.11543936 
## sample estimates: 
##     prop 1     prop 2  
## 0.20211757 0.09564045 

このように3群いずれにいても有意差が見られました。しかし、今回の場合は先程見た若者とは逆の有意差で今回は実際の人口比に比べて50代以上の感染者比が小さいことがわかります。（繰り返しにはなりますが、これだけで判断できるものではありません。）

SIRモデルを用いたシミュレーション

次にSIRモデルと呼ばれる感染者数の予測の際によく用いられるモデルを用いて感染者数の予測を行います。
SIRモデルについてはこちらのページをご覧ください。

これらの微分方程式から$\frac{dI}{dt}=I(\beta-\gamma)>0$の時に感染は拡大します。
従って、$R_0=\frac{\beta S}{\gamma}>1$の時に感染は拡大するといえます。この$R_0$が基本再生産数と呼ばれ、1人の患者が平均何名に新たに罹患させるかを表す値です。つまり、終息のためにはこの$R_0$をできるだけ小さくすることが求められます。

以下ではこの$R_0$を主役として考えていきます。

まず、なんとかして今あるデータから$R_0$を導出したいと思います。$\frac{dI}{dt}=I(\beta S-\gamma)$より患者の数は$(\beta S-\gamma)$の等比級数で増大します。ここで簡単のために$S \simeq 1$とします。(現段階では多くの人たちが免疫を獲得していないので割と妥当な仮定です。)
すると感染者数は$(\beta - \gamma)$の等比級数で増大します。
以上のことをまとめると、発症者数を$n$,確定日を$t$としたときに$log(n)=at+b$と線形回帰した時の$a$が$(\beta - \gamma)$となるということです。(患者の数は$(\beta S-\gamma)$の等比級数で増大するので)

実際に線形回帰を行ってみます。

n<-lm(log(通し)~確定日,data=cases)
plot(cases$確定日,log(cases$通し),xlab="確定日",ylab="log(感染者数)")
n
abline(n)

こうすると係数が以下のように出てきます。

##  ## Call: 
## lm(formula = log(通し) ~ 確定日, data = cases) 
##  ## Coefficients: 
## (Intercept)       確定日   
##  -1.480e+03    8.109e-02 

また、回帰を行った時のグラフは以下のようになります。

以上のことをまとめると$\beta-\gamma \simeq 0.08$となります。
ここで$\gamma$について考えます。今後の解析を簡単にするために$\gamma \simeq 0.1$とします。つまり、新型コロナウイルスに罹患してから回復まで約10日間かかるとしました。(ここの仮定は結構乱暴かもしれません…)
以上の仮定から$\beta \simeq 0.18$となります。
これらを用いると$S$を用いて基本再生産数$R_0$は$1.8S$と表すことが出来ます。($R_0=\frac{\beta S}{\gamma}$から)

これでSIRモデルを用いてシミュレーションをする準備が整いました。
早速、以上の情報を用いてSIRモデルを用いたシミュレーションを行います。
以下の解き方はこちらを参考にさせていただきました。

library(dplyr)
library(tidyr)
library(ggplot2)
library(gridExtra)
#国内ではじめて感染者が確認されてから1年間の感染者数の推移を予測する。
sir<-function(beta, gamma){
  a<-numeric(366)
  b<-numeric(366)
  c<-numeric(366)
  b[1]<-1e-6
  c[1]<-0
  a[1]<-1-b[1]
  for (i in 1:365) {
    a[i+1]<-a[i]-beta*a[i]*b[i]
    b[i+1]<-b[i]+beta*a[i]*b[i]-gamma*b[i]
    c[i+1]<-c[i]+gamma*b[i]
  }
  data.frame(days=0:365, s=a, i=b, r=c)
}
theme_set(theme_classic(base_size = 18, base_family = "Helvetica"))
result1<-sir(0.1809,0.1)# beta=0.1809,gamma=0.1で上の関数を動かす。
pivot_longer(result1, cols=-days,
             names_to="situation", values_to="ratio") %>%
  ggplot(aes(days,ratio, color=situation, linetype=situation)) +
  geom_path() +
  theme_classic()+
  theme_bw()

すると、以下のような図が得られます。

このように外出制限などなにも対策を行わなかった場合、最も流行しているときには日本国民の12.5%が罹患していることになります。日本の人口を1憶人として計算するとこれは実に12500000人が感染していることになります。(執筆時点で感染者数の増加は抑えらえつつあるのでこのような事態は避けられたと思います。)

次に厚生労働省の専門家チームが発表している「このまま何の対策しなければ約40万人の人が亡くなる」という言説について妥当性を確かめたいと思います。そのために次の仮定を置きたいと思います。

• 罹患者数のうち0.5%の方が亡くなるとする。(罹患者数の2%の方が重症化するということより)
• 集団免疫を獲得するころには75%の人が感染する。(これは先ほど示した図からの判断です。)

以上から簡単に次のような計算が出来ます。$100000000 \times 0.75 \times 0.005 =375000$となる。さらに、最流行時には病床の数も不足することにより、もっと多くの方が亡くなると考えると約40万人の方が亡くなると予想するのは妥当であることがわかりました。(もちろん、今は外出自粛などの効果で抑えられてはいます。)

今までのものは全て外出自粛などの対策をとらなかった場合のものです。それでは、対策をした場合はどうなるのでしょうか。これについてシミュレーションをしてみましょう。
以下のシュミレーションでは緊急事態宣言が発出された4月7日という日本国内においてはじめて感染者が確認されて84日目までを$R_0=1.8$と計算して、それ以降の$R_0$をいろいろいじってみたいと思います。
$R_0$のいじり方ですが今回は以下のようにいじります。

• $R_0=\frac{\beta S}{\gamma}$なので、$R_0$を減らすためには$\gamma$を増やすか$\beta$を減らすことになります。($S$は未感染者の割合を表すものであり、今考えているパラメタは$\gamma,\beta$です。)
• しかし、$\gamma$を増やすということは回復率を上げることを意味します。このためには病床の数を増やしたり、ワクチンの開発が必要です。これらは短期的に達成することはできないので、今回は$R_0$をいじるということは$\beta$をいじることを意味するという下でシミュレーションを行います。

さきに、先ほどのシミュレーションで用いたSIRモデルでは詳細な時間の範囲を指定できなかったので、時間の範囲を指定できるように先ほどの関数を書き換えます。
そして、この新たに定義した関数を用いて日本国内初の感染者確認から84日目までを$R_0=1.8$の下で計算します。これはさきほどの何も対策をしなかった場合のシミュレーションの84日目までと全く同じです。

sir2<-function(beta,gamma,MAXT,b0,c0){
  a0<-1-b0-c0
  a<-numeric(1+MAXT)
  b<-numeric(1+MAXT)
  c<-numeric(1+MAXT)
  a[1]<-a0
  b[1]<-b0
  c[1]<-c0
  for (i in 1:MAXT) {
    a[i+1]<-a[i]-beta*a[i]*b[i]
    b[i+1]<-b[i]+beta*a[i]*b[i]-gamma*b[i]
    c[i+1]<-c[i]+gamma*b[i]
  }
  data.frame(t=0:MAXT,s=a,i=b,r=c)
}
simdata<-sir2(0.1809,0.1,84,1e-6,0)
pivot_longer(simdata, cols=-t,
             names_to="type", values_to="fraction") %>%
  ggplot(aes(t, fraction, color=type, linetype=type)) +
  geom_path() +
  ylim(0,0.0025)+
  theme_bw()

次に実際に85日目以降の$\beta$を変更してみます。($R_0$をいじることと同値)

#85日目以降の\beta=0.15としたとき
simdata<-sir2(0.15,0.1,241,6.864532e-04,8.4868e-04)
pivot_longer(simdata, cols=-t,
             names_to="type", values_to="fraction") %>%
  ggplot(aes(t, fraction, color=type, linetype=type)) +
  geom_path() +
  ylim(0,0.25)+
  theme_bw()

#85日目以降の\beta=0.10としたとき
simdata<-sir2(0.1,0.1,241,6.864532e-04,8.4868e-04)
pivot_longer(simdata, cols=-t,
             names_to="type", values_to="fraction") %>%
  ggplot(aes(t, fraction, color=type, linetype=type)) +
  geom_path() +
  ylim(0,0.005)+
  theme_bw()

#85日目以降の\beta=0.08としたとき
simdata<-sir2(0.08,0.1,241,6.864532e-04,8.4868e-04)
pivot_longer(simdata, cols=-t,
             names_to="type", values_to="fraction") %>%
  ggplot(aes(t, fraction, color=type, linetype=type)) +
  geom_path() +
  ylim(0,0.001)+
  theme_bw()

以上のように複数の$\beta$についてシミュレーションを行いました。($t$軸は緊急事態宣言発出日より数えた日数)
以上のグラフを見たらわかるように$\beta\leq0.1$にすることで緊急事態宣言発出時の感染者を維持、さらには減少させることが出来ます。(これは、$\gamma\simeq0.1$としていることよりもわかります。$S\leq1$より$\beta\leq0.1$において必ず$R_0\leq1$となります。)それ以上の値ではピーク時の感染者は何も対策を行い時に比べては減少はするが、十分でないです。

それでは、$\beta$を実際に$0.18$から$0.1$にするためにはどうすればよいのでしょう。$\beta$の定義に立ち戻るとこれは$S,I$が接触した時に$S$が感染する割合だと解釈でき、その低下のためにはそもそもの接触数を減らすことが第一であることがわかります。それでは、具体的にどれほど減少させればよいのでしょう。$\beta$を決定付けているのは様々な要因が考えられるため、今回のシミュレーションの結果のみから判断はできないですが、$0.1/0.18\simeq0.556$と算出できます。つまり、$\beta$の変化を人と人の接触のみに依存するものだとすると、約$45$%接触を減らすと、所望の$\beta$が得られることになります。しかし、これは何度も言うように人と人の接触しか考えていない上に、$\beta$の変化には他の要因もあるのでこれは正しいとは言えません。さらには、今回は$\gamma=0.1$の定数として扱っていて、これも誤差を生み出す原因となっています。実際に、報道で伝えられている通り推奨されているのは接触$8$割減です。

参考文献

https://jag-japan.com/covid19map-readme/
https://www.stat.go.jp/index.html
http://ranalytics.blog.fc2.com/blog-entry-32.html
https://ja.wikipedia.org/wiki/SIR%E3%83%A2%E3%83%87%E3%83%AB
https://qiita.com/kota9/items/c08373d4ef8230922aa7