◆【新型コロナ】ECDCのデータから国別の実効再生産数を計算するRのコード例です

ECDCが公開している新型コロナウイルスの感染確認者数のデータを使って、各国の実効再生産数を計算するコードの例です。

日付のデータを付加するために、コードが長くなっていますが、日付データを付加しなければもっとシンプルにすることができます。

古いノートパソコンを使っていますが、約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")

          ------------------------------------------------------------------------------

 

------------------------------------------------------------------------------

-------------------------------------------------------------------

--------------------------------------------------------------------