流行予測で用いているSIRforecastのカスタムスクリプト。 行に日付、列に地域、値に人口比の累計感染者数をまとめたデータフレームを受け取ると、recoverdays（デフォルトは14日）でこれら感染者が回復する前提で、averagedays（デフォルトは14日）分の前方移動平均によるβ（感染率）、γ（回復率）を算定、SIRモデルに基づくpredictdays（デフォルトは360日）分の予測値を、直近、1週間前、2週間前の3つの時点での感染者数の伸び率に基づくβ及びγから算出。
SIRforecast<- function(df=df,
recoverdays=14,
averagedays=14,
para_average=14,
predictdays=360){
#日付と分離
df$date<- as.Date(df$date)
df<- df[1:(nrow(df)-1),]
date<- df[,1]
para_date<- df$date[(para_average):length(df$date)]
#S,I,Rの算出
work<- df[,2:ncol(df)]
S<- 1-work
dS<- data.frame(diff(as.matrix(S)))
dS[2:(nrow(dS)+1),]<- dS
dS[1,]<- NA
R<- data.frame(
apply(work,2,function(x){lag(x,n=recoverdays)}))
dR<- data.frame(diff(as.matrix(R)))
dR[2:(nrow(dR)+1),]<- dR
dR[1,]<- NA
I<- 1-S-R
dI<- data.frame(diff(as.matrix(I)))
dI[2:(nrow(dI)+1),]<- dI
dI[1,]<- NA
cumI<- work
dcumI<- data.frame(diff(as.matrix(cumI)))
dcumI[2:(nrow(dcumI)+1),]<- dcumI
dcumI[1,]<- NA
#β、γの算出、移動平均化
for_para_I<- apply(I,2,function(x){ifelse(x==0,NA,x)})
# raw_beta<- -dS/(apply(S,2,function(x){lag(x)})*
# apply(for_para_I,2,function(x){lag(x)}))
raw_beta<- -dS/(S*for_para_I)
ave_beta<- data.frame(
apply(raw_beta,2,
function(x){roll_mean(x,n=averagedays,align="right")}))
ave2_beta<- data.frame(
apply(raw_beta,2,
function(x){roll_mean(x,n=para_average,align="right")}))
beta<- S
beta[1:(averagedays-1),]<- NA
beta[averagedays:nrow(beta),]<- ave_beta
# raw_gamma<- dR/apply(for_para_I,2,function(x){lag(x)})
raw_gamma<- dR/apply(for_para_I,2,function(x){lag(x)})
ave_gamma<- data.frame(
apply(raw_gamma,2,
function(x){roll_mean(x,n=averagedays,align="right")}))
ave2_gamma<- data.frame(
apply(raw_gamma,2,
function(x){roll_mean(x,n=para_average,align="right")}))
gamma<- S
gamma[1:(averagedays-1),]<- NA
gamma[averagedays:nrow(gamma),]<- ave_gamma
beta_recent<- data.frame(
t(apply(beta,2,function(x){last(x)})))
beta_2weeksago<- data.frame(
t(apply(beta,2,function(x){nth(x,n=nrow(beta)-14)})))
beta_1weeksago<- data.frame(
t(apply(beta,2,function(x){nth(x,n=nrow(beta)-7)})))
gamma_recent<- data.frame(
t(apply(gamma,2,function(x){last(x)})))
gamma_2weeksago<- data.frame(
t(apply(gamma,2,function(x){nth(x,n=nrow(gamma)-14)})))
gamma_1weeksago<- data.frame(
t(apply(gamma,2,function(x){nth(x,n=nrow(gamma)-7)})))
#変数出力
out_beta<- transform(ave2_beta,date=para_date)
out_beta<- out_beta[,c(ncol(out_beta),1:(ncol(out_beta)-1))]
out_gamma<- transform(ave2_gamma,date=para_date)
out_gamma<- out_gamma[,c(ncol(out_gamma),1:(ncol(out_gamma)-1))]
ave2_Rt<- ave2_beta/ave2_gamma
out_Rt<- transform(ave2_Rt,date=para_date)
out_Rt<- out_Rt[,c(ncol(out_Rt),1:(ncol(out_Rt)-1))]
#初期値の設定
init_S_recent<- data.frame(apply(S,2,function(x){last(x)}))
init_S_2weeksago<- data.frame(
apply(S,2,function(x){nth(x,n=nrow(S)-14)}))
init_S_1weeksago<- data.frame(
apply(S,2,function(x){nth(x,nrow(S)-7)}))
init_I_recent<- data.frame(apply(I,2,function(x){last(x)}))
init_I_2weeksago<- data.frame(
apply(I,2,function(x){nth(x,n=nrow(I)-14)}))
init_I_1weeksago<- data.frame(
apply(I,2,function(x){nth(x,nrow(I)-7)}))
init_newI_recent<- data.frame(apply(dcumI,2,function(x){last(x)}))
init_newI_2weeksago<- data.frame(
apply(dcumI,2,function(x){nth(x,n=nrow(dcumI)-14)}))
init_newI_1weeksago<- data.frame(
apply(dcumI,2,function(x){nth(x,nrow(dcumI)-7)}))
init_R_recent<- data.frame(apply(R,2,function(x){last(x)}))
init_R_2weeksago<- data.frame(
apply(R,2,function(x){nth(x,n=nrow(R)-14)}))
init_R_1weeksago<- data.frame(
apply(R,2,function(x){nth(x,nrow(R)-7)}))
#結果格納データフレーム
pred_I_recent<- data.frame(
matrix(NA,nrow=predictdays,ncol=ncol(S)))
colnames(pred_I_recent)<- colnames(S)
pred_I_2weeksago<- pred_I_recent
pred_I_1weeksago<- pred_I_recent
pred_newI_recent<- pred_I_recent
pred_newI_2weeksago<- pred_I_recent
pred_newI_1weeksago<- pred_I_recent
pred_R_recent<- pred_I_recent
pred_R_2weeksago<- pred_I_recent
pred_R_1weeksago<- pred_I_recent
pred_S_recent<- pred_I_recent
pred_S_2weeksago<- pred_I_recent
pred_S_1weeksago<- pred_I_recent
pred_I_recent[1,]<- matrix(
t(init_I_recent),nrow=1,ncol=ncol(S))
pred_I_2weeksago[1,]<- matrix(
t(init_I_2weeksago),nrow=1,ncol=ncol(S))
pred_I_1weeksago[1,]<- matrix(
t(init_I_1weeksago),nrow=1,ncol=ncol(S))
pred_newI_recent[1,]<- matrix(
t(init_newI_recent),nrow=1,ncol=ncol(S))
pred_newI_2weeksago[1,]<- matrix(
t(init_newI_2weeksago),nrow=1,ncol=ncol(S))
pred_newI_1weeksago[1,]<- matrix(
t(init_newI_1weeksago),nrow=1,ncol=ncol(S))
pred_R_recent[1,]<- matrix(t(init_R_recent),nrow=1,ncol=ncol(S))
pred_R_2weeksago[1,]<- matrix(
t(init_R_2weeksago),nrow=1,ncol=ncol(S))
pred_R_1weeksago[1,]<- matrix(
t(init_R_1weeksago),nrow=1,ncol=ncol(S))
pred_S_recent[1,]<- matrix(t(init_S_recent),nrow=1,ncol=ncol(S))
pred_S_2weeksago[1,]<- matrix(
t(init_S_2weeksago),nrow=1,ncol=ncol(S))
pred_S_1weeksago[1,]<- matrix(
t(init_S_1weeksago),nrow=1,ncol=ncol(S))
#予測
for(i in 2:predictdays){
pred_S_recent[i,]<- pred_S_recent[(i-1),]-
pred_S_recent[(i-1),]*pred_I_recent[(i-1),]*beta_recent
pred_I_recent[i,]<- pred_I_recent[(i-1),]+
pred_S_recent[(i-1),]*pred_I_recent[(i-1),]*beta_recent-
pred_I_recent[(i-1),]*gamma_recent
pred_R_recent[i,]<- pred_R_recent[(i-1),]+
pred_I_recent[(i-1),]*gamma_recent
pred_newI_recent[i,]<-
pred_S_recent[(i-1),]*pred_I_recent[(i-1),]*beta_recent
pred_S_2weeksago[i,]<- pred_S_2weeksago[(i-1),]-
pred_S_2weeksago[(i-1),]*
pred_I_2weeksago[(i-1),]*beta_2weeksago
pred_I_2weeksago[i,]<- pred_I_2weeksago[(i-1),]+
pred_S_2weeksago[(i-1),]*
pred_I_2weeksago[(i-1),]*beta_2weeksago-
pred_I_2weeksago[(i-1),]*gamma_2weeksago
pred_R_2weeksago[i,]<- pred_R_2weeksago[(i-1),]+
pred_I_2weeksago[(i-1),]*gamma_2weeksago
pred_newI_2weeksago[i,]<-
pred_S_2weeksago[(i-1),]*pred_I_2weeksago[(i-1),]*
beta_2weeksago
pred_S_1weeksago[i,]<- pred_S_1weeksago[(i-1),]-
pred_S_1weeksago[(i-1),]*
pred_I_1weeksago[(i-1),]*beta_1weeksago
pred_I_1weeksago[i,]<- pred_I_1weeksago[(i-1),]+
pred_S_1weeksago[(i-1),]*
pred_I_1weeksago[(i-1),]*beta_1weeksago-
pred_I_1weeksago[(i-1),]*gamma_1weeksago
pred_R_1weeksago[i,]<- pred_R_1weeksago[(i-1),]+
pred_I_1weeksago[(i-1),]*gamma_1weeksago
pred_newI_1weeksago[i,]<-
pred_S_1weeksago[(i-1),]*pred_I_1weeksago[(i-1),]*
beta_1weeksago
}
#日付列作成
date_recent<- seq(last(df$date),
by="1 day",length.out=predictdays)
date_2weeksago<- seq(last(df$date)-days(14),
by="1 days",length.out=predictdays)
date_1weeksago<- seq(last(df$date)-days(7),
by="1 days",length.out=predictdays)
#日付の付与
pred_S_recent<- transform(pred_S_recent,date=date_recent)
pred_I_recent<- transform(pred_I_recent,date=date_recent)
pred_newI_recent<- transform(pred_newI_recent,date=date_recent)
pred_R_recent<- transform(pred_R_recent,date=date_recent)
pred_S_2weeksago<- transform(
pred_S_2weeksago,date=date_2weeksago)
pred_I_2weeksago<- transform(
pred_I_2weeksago,date=date_2weeksago)
pred_newI_2weeksago<- transform(
pred_newI_2weeksago,date=date_2weeksago)
pred_R_2weeksago<- transform(
pred_R_2weeksago,date=date_2weeksago)
pred_S_1weeksago<- transform(
pred_S_1weeksago,date=date_1weeksago)
pred_I_1weeksago<- transform(
pred_I_1weeksago,date=date_1weeksago)
pred_newI_1weeksago<- transform(
pred_newI_1weeksago,date=date_1weeksago)
pred_R_1weeksago<- transform(
pred_R_1weeksago,date=date_1weeksago)
#ワイドからロングへ
dfoutpara01<- gather(out_beta,key="region",value="value",-date)
dfoutpara01<- transform(dfoutpara01,param="beta")
dfoutpara01<- transform(dfoutpara01,type=NA)
dfoutpara02<- gather(out_gamma,key="region",value="value",-date)
dfoutpara02<- transform(dfoutpara02,param="gamma")
dfoutpara02<- transform(dfoutpara02,type=NA)
dfoutpara03<- gather(out_Rt,key="region",value="value",-date)
dfoutpara03<- transform(dfoutpara03,param="Rt")
dfoutpara03<- transform(dfoutpara03,type=NA)
dfout01<- gather(pred_S_recent,key="region",value="value",-date)
dfout01<- transform(dfout01,param="S")
dfout01<- transform(dfout01,type="recent")
dfout02<- gather(pred_I_recent,key="region",value="value",-date)
dfout02<- transform(dfout02,param="I")
dfout02<- transform(dfout02,type="recent")
dfout04<- gather(pred_newI_recent,key="region",value="value",-date)
dfout04<- transform(dfout04,param="newI")
dfout04<- transform(dfout04,type="recent")
dfout03<- gather(pred_R_recent,key="region",value="value",-date)
dfout03<- transform(dfout03,param="R")
dfout03<- transform(dfout03,type="recent")
dfout11<- gather(pred_S_2weeksago,key="region",value="value",-date)
dfout11<- transform(dfout11,param="S")
dfout11<- transform(dfout11,type="2weeks")
dfout12<- gather(pred_I_2weeksago,key="region",value="value",-date)
dfout12<- transform(dfout12,param="I")
dfout12<- transform(dfout12,type="2weeks")
dfout14<- gather(pred_newI_2weeksago,key="region",value="value",-date)
dfout14<- transform(dfout14,param="newI")
dfout14<- transform(dfout14,type="2weeks")
dfout13<- gather(pred_R_2weeksago,key="region",value="value",-date)
dfout13<- transform(dfout13,param="R")
dfout13<- transform(dfout13,type="2weeks")
dfout21<- gather(pred_S_1weeksago,key="region",value="value",-date)
dfout21<- transform(dfout21,param="S")
dfout21<- transform(dfout21,type="1weeks")
dfout22<- gather(pred_I_1weeksago,key="region",value="value",-date)
dfout22<- transform(dfout22,param="I")
dfout22<- transform(dfout22,type="1weeks")
dfout24<- gather(pred_newI_1weeksago,key="region",value="value",-date)
dfout24<- transform(dfout24,param="newI")
dfout24<- transform(dfout24,type="1weeks")
dfout23<- gather(pred_R_1weeksago,key="region",value="value",-date)
dfout23<- transform(dfout23,param="R")
dfout23<- transform(dfout23,type="1weeks")
#結合・出力
dfout<- rbind(dfout01,dfout02,dfout03,dfout04,
dfout11,dfout12,dfout13,dfout14,
dfout21,dfout22,dfout23,dfout24,
dfoutpara01,dfoutpara02,dfoutpara03)
dfout
}