日付のデータを付加するために、コードが長くなっていますが、日付データを付加しなければもっとシンプルにすることができます。
古いノートパソコンを使っていますが、約200カ国の実効再生産数を十数秒程度で計算することができます。
計算で得られたデータは、ダッシュボードで利用しています。
なお、実効再生産数の計算には、「Improved inference of time-varying reproduction numbers during infectious disease outbreaks」という論文で紹介されているRのパッケージ「EpiEstim」を利用しています。
なお、発症間隔のパラメータは、平均を4.8、標準偏差を2.3に設定しています。
【Rコードの例】
library(EpiEstim)
df_ECDC <- read.csv("https://opendata.ecdc.europa.eu/covid19/casedistribution/csv", na.strings = "", fileEncoding = "UTF-8-BOM",stringsAsFactors = FALSE)
geo_list <- unique(df_ECDC$countriesAndTerritories)
df_ECDCtemp <- df_ECDC
df_ECDCtemp1 <- NULL
df_ECDCtemp2 <- NULL
temp_R <- NULL
temp_Rt <- NULL
temp_Date <- NULL
temp_Case <- NULL
temp_DC <- NULL
df_DC <- NULL
df_Rt <- NULL
temp_notcal <- NULL
df_notcal <- NULL
###繰り返しの処理で、多くの国の計算を行います。
for (i in seq_along(geo_list))
{
df_ECDCtemp1 <- subset(df_ECDCtemp,df_ECDCtemp$countriesAndTerritories==geo_list[i])
df_ECDCtemp2 <- subset(df_ECDCtemp1,df_ECDCtemp1$cases >= 0)
if (length(df_ECDCtemp2$cases) >= 70)
{rt_parametric_si <- estimate_R(df_ECDCtemp2$cases,method = "parametric_si",config = make_config(list(mean_si = 4.8,
std_si = 2.3)))
temp_R <- round(rt_parametric_si$R,2)
temp_Rt <- mutate(temp_R,countriesAndTerritories=geo_list[i])
Numa <- nrow(rt_parametric_si$R)
temp_date <- as.data.frame(df_ECDCtemp1$Date)
Numb <- nrow(temp_date)
Numc <- Numb-Numa
temp_date <- temp_date[-c(1:Numc),]
temp_Rt <- mutate(temp_Rt,Days=seq(from=1,to=nrow(rt_parametric_si$R), by=1))
temp_Rt <- mutate(temp_Rt,Date=temp_date)
temp_Date1 <- matrix(rt_parametric_si$dates,ncol=1)
colnames(temp_Date1) <- "Days"
temp_Case <- matrix(rt_parametric_si$I,ncol=1)
colnames(temp_Case) <- "Cases"
temp_DC <- cbind(temp_Date1,temp_Case)
temp_DC <- as.data.frame(temp_DC)
temp_DC <- mutate(temp_DC,countriesAndTerritories=geo_list[i])
df_DC <- rbind(df_DC,temp_DC)
df_Rt <- rbind(df_Rt,temp_Rt)}
else {temp_notcal <- as.data.frame(geo_list[i])
temp_notcal <- mutate(temp_notcal,Under70days=nrow(df_ECDCtemp1))
colnames(temp_notcal) <- c("countriesAndTerritories","Under70days")
df_notcal <- rbind(df_notcal,temp_notcal)
}
}
###グラフで「１」の線を引くためのデータを追加します
Const <- c(rep(1,nrow(df_Rt)))
Const <- as.data.frame(Const)
colnames(Const) <- c("C")
df_Rt <- cbind(Const,df_Rt)
write.csv(df_Rt,"Rt.csv",fileEncoding = "UTF8")
------------------------------------------------------------------------------
------------------------------------------------------------------------------
--------------------------------------------------------------------
