感染症モデル入門（4)：統計言語[R]によるSIRモデルのシミュレーション

１．SIRモデルを統計言語[R]でシミュレーションする。

以下にRのプログラムを示す。なるべくシンプルなコーディングを目指したので、統計言語[R]の初歩的な知識があれば、コメント(＃)行をみれば、何をやっているのかがわかるであろう。
このプログラムは、自分のPCに[R]が既にインストールしてあれば、その上で動く。あるいは、手軽に試してみたい場合には、ネット上にオンラインで[R]を走らせるサイトがいくつかある。たとえば、

1)　ここcompile R onlineをクリックして、以下のプログラムをコピペして、左下にある[run it]をクリックするか、自分のPCでF8を押せばよい。あるいは、
2)　ここpaiza.ioでも同様なことができる。

試してみてはいかがだろうか？

# ==================================================
#  [R]によるSIRモデルの計算 ＃#==================================================
#####　シミュレーション期間の設定とS,I,Rベクトルの定義
No_of_Days=50; N1=No_of_Days+1; S=I=R=R_t=numeric(N1); 
##### S,I,Rの初期値とパラメータ値の設定
S[1]=100; I[1]=1; R[1]=0; beta=0.01; gamma=0.1
#####　SIRモデルの繰り返し計算
for (t in 1:No_of_Days){
    S[t+1] = S[t] - beta*S[t]*I[t]
    I[t+1] = I[t] + beta*S[t]*I[t] - gamma*I[t]
    R[t+1] = R[t] + gamma*I[t]
    R_t[t+1] = S[t+1]*beta/gamma
}
##### S,I,Rのグラフを描く（S,I,Rと横軸に日数をひとまとめのデータフレームとしてから）
output = data.frame( S, I, R, date=1:N1 )
library(ggplot2)
ggplot( output, aes(date)) +
  geom_line(aes(y = S, colour = "S")) +
  geom_line(aes(y = I, colour = "I")) +
  geom_line(aes(y = R, colour = "R"))
##### 再生産数 R_tのグラフを描く
plot(R_t, type="l", lwd=2, xaxt="n")
# ==================================================

2.  S,I,R の推移をグラフでしめす

上のプログラムを実行すれば、次のようSIRと実行再生産数R_tの推移をえることができる。
自分なりのシナリオを描いて、政策シミュレーションをやってみたい場合には、上のプログラムの

##### S,I,Rの初期値とパラメータ値の設定

S[1]=100; I[1]=1; R[1]=0; beta=0.01; gamma=0.1

で、それぞれの初期値S[1]=100; I[1]=1; R[1]=0;　と　感染率と回復率のパラメータを変化させてることで可能である。

統計言語[R]によるSIRモデルの計算結果

Rによる実行再生算数：R_t