新型コロナウイルスの感染者数を予測する で感染者数のトレンドと2週間後の予測を出してみました。
この結果を利用して、実行再生産数を計算してみます。
実行再生産数とは、一人の感染者が何人に感染させたかを表す数値です。
１以上なら感染は拡大していき、1未満ならいずれ収束していくことになります。
実行再生産数を求める方法は、下記を使わせていただきました。
※上記サイトでは発症から検査確定までの期間（感染期間）に感染していると想定したモデルを使用していますが、潜伏期間中も感染はあると考え、その点のみ変更して利用しています。
それでは東京都の結果です。
実効再生産数（緑線）は感染者数そのものではなく、予測値（青線）の値から算出しました。
感染者数から直接算出すると、がたつきの多いグラフになり、傾向が読み取りづらかったためです。（下図のようになります）
実行再生産数は、13日（感染してから検査で確定するまでの期間）先までの感染者数から遡って割り出しています。そのため、現在から13日前から徐々に予測感染者数に基づいた値になっていくため、トレンドに相関した曲線になっていきます。
まず、感染者数についてですが、今日（4/2）の97人が効いて、感染者数の予測は悪化しています。下記は前回（3/31）に作成した図。
4/10の予測数は、前回は約110人、今回は約130人となっています。
感染者数は増加傾向ですが、実効再生産数は１を若干上回る程度であり、自粛等の効果により、ある程度抑制されていると思います。
しかし、実効再生産数が１を下回らない限り、収束に向かうことはありません。
たとえ実効再生産数＝１だったとしても、感染者は一定のペースで増え続けることになります。同じペースで回復していく人がいれば均衡は保たれるわけですが、感染してから入院するまでの期間より、入院から退院までの方が時間がかかるわけですから、入院患者は増える一方です。
専門家会議の先生が、「オーバーシュートの前に医療崩壊が起こる」と言ったのは、こういうことかな、と思います。
次に全国。
こちらは東京より希望が持てる感じですが、東京ほど地方では増えていないということであり、やはり各都市部や県で見た方が実態を捉えられるのかもしれません。（次回はそうしよう）
(4/5修正)
4/2（この記事作成当日）の件数が少なかったため、正しいグラフになっていませんでした。
当日のデータは集計結果が徐々に反映されるようなので、注意が必要ですね。
修正したグラフでは、東京と類似した結果となっています。
library(dplyr)
library(KFAS)
setwd("E:/F/Statistics/Corona/Src")
#感染者数推定：２変量平滑化トレンドモデル
areaname = "東京都"
inf_cnt = 13 #感染してから確定するまでの期間（潜伏期間＋発症期間）
predict_cnt = 30 #予測期間
y_scale=c(0,200)
if (areaname=="全国") y_scale=c(0,300)
rawdata <- read.csv("COVID-19.csv", fileEncoding = "UTF-8-BOM")
areadata <- rawdata %>% filter(受診都道府県==areaname)
if (nrow(areadata) == 0) areadata = rawdata
cntdata <- areadata %>% group_by(確定日YYYYMMDD, .drop=FALSE) %>% summarise(count = n())
colnames(cntdata) <- c("date","count")
cntdata[[1]] <- as.Date(cntdata[[1]])
cntdata <- cntdata %>% filter(date > "2020-02-06")
cntdata <- cntdata[order(cntdata$date),]
#最新値が０ならNaNにする
rcnt = length(cntdata$count)
if (cntdata$count[rcnt] == 0) cntdata$count[rcnt] = NA
#値のない日を挿入
dts = seq(min(cntdata$date), max(cntdata$date) + predict_cnt, "day")
dtframe = data.frame(date=dts)
cntdata <- left_join(dtframe, cntdata ,by="date")
Pobj <- ts(cntdata$count)
mod <- SSModel(Pobj ~ SSMtrend(2, Q = list(0, NA)), H = NA)
#fit_method = "SANN"
fit_method = "BFGS"
fit <- fitSSM(mod, numeric(2), method=fit_method)
kfs <- KFS(fit$model)
conf <- predict(fit$model, interval="confidence", level=0.95)
pre <- predict(fit$model, interval="prediction", level=0.95)
logLik <- kfs$logLik - sum(kfs$Finf>0) * log(2*pi)/2
AIC <- -2*logLik + 2*( 2 + 2 )
#図の描画
par(oma = c(0, 0, 0, 2))
par(xaxt="n")
e = length(Pobj) - predict_cnt/2
showdts = dts[1:e]
#plot(showdts,Pobj[1:e],type="l",col="red",lwd=2,xlab="",ylab=paste("１日当たり新規感染者数(",areaname,")"),ylim=c(0,200))
barplot(Pobj[1:e],col="orange",xlab="",ylab=paste("１日当たり新規感染者数(",areaname,")"),ylim=y_scale,axes=FALSE)
par(new=T)
plot(showdts,pre[,1][1:e],type="l",col="blue",lwd=2,xlab="",ylab="",ylim=y_scale, yaxs="i")
#lines(showdts,pre[,1][1:e],col="blue",lwd=2)
s = length(Pobj) - predict_cnt
lines(showdts[s:e],pre[s:e,2],col="blue",lty=2)
lines(showdts[s:e],pre[s:e,3],col="blue",lty=2)
par(xaxt = "s")
axis.Date(1, at=seq(min(showdts), max(showdts), "7 day"), format="%m/%d", tck=1.0,lty="dotted")
axis(side=2,tck=1.0,lty="dotted")
#実効再生産数
ee = length(Pobj)
cntdata$count[(s+1):ee] = pre[(s+1):ee,1] #実値と予測値を結合
R0obj = lapply(1:(ee-inf_cnt-1), function(i) {
inf_sum = sum(pre[(i+1):(i+1+inf_cnt),1], na.rm=TRUE)
if (inf_sum == 0) return(NA)
return(inf_cnt * pre[i+inf_cnt,1] / inf_sum)
}
)
#図の追加描画
par(new=T)
plot(showdts,R0obj[1:e],type="l",col="forestgreen",lwd=2,xlab="",ylab="",ylim=c(0,4), yaxs="i",axes=FALSE)
abline(h=1, col='red', lwd=1, lty=2)
mtext("実効再生産数",side = 4, line = 3)
axis(4)
legend("topleft",legend = c("新規感染者数","トレンド（予測）","実効再生産数（予測）"),lwd=3,col=c("orange","blue","forestgreen"))
「新型コロナウイルス感染者数マップ」
