Calculate Estimated R in Japan
In this article, I will show how to plot an interactive effective reproduction number line plot using R and {echarts4r}. However, since I do not know any epidemic theory, I will not explain too much about the model.
The same code also been used in COVID-19 BULLETIN BOARD.
1. Install Packages
First, install all the packages we need. Some minor bugs in {echarts4r} has been fixed in development version but has not been release to CRAN, so we need to install the Github version here.
2. Create the Dataset for Plotting
# Load packages
library(echarts4r)
library(incidence)
library(EpiEstim)
library(data.table)
# Read the dataset from `2019-ncov-japan`
byDate <-
fread(
"https://raw.githubusercontent.com/swsoyee/2019-ncov-japan/master/50_Data/byDate.csv"
)
byDate[is.na(byDate)] <- 0
# Paramsters refer to https://www.ijidonline.com/article/S1201-9712(20)30119-3/fulltext
mean_si <- 4.6
std_si <- 2.6
# Store all names of prefecture
pref <- colnames(byDate)[2:48]
# For storing all echarts4r object in to one htmlwidgets objects
output_list <- list()
#' Function for generate a dataset contains dates and incidences, and Rt
#'
#' @param dataset dataset from `2019-ncov-japan`.
#' @param pref prefecture name.
#'
#' @return a `data.table` for ploting.
createEstimateFromSource <- function(dataset, pref) {
incid <-
as.incidence(
rowSums(dataset[, pref, with = FALSE]),
dates = as.Date(as.character(dataset$date), format = "%Y%m%d")
)
res <- suppressMessages(suppressWarnings(
estimate_R(incid,
method = "parametric_si",
config = make_config(list(
mean_si = mean_si,
std_si = std_si,
t_end = max(incid$dates)
))
)
))
dt <- as.data.table(res$R)
cols <- colnames(dt)
dt[, (cols) := lapply(.SD, function(x) {
return(round(x, 2))
}), .SDcols = cols]
dt$dates <- res$dates[res$R$t_end]
dt$Incidence <- res$I[res$R$t_end]
dt
}
# Test
tokyo <- createEstimateFromSource(byDate, "東京")
tokyo## t_start t_end Mean(R) Std(R) Quantile.0.025(R) Quantile.0.05(R)
## 1: 2 8 1.56 0.90 0.32 0.43
## 2: 3 9 0.93 0.66 0.11 0.16
## 3: 4 10 0.90 0.64 0.11 0.16
## 4: 5 11 0.97 0.68 0.12 0.17
## 5: 6 12 1.07 0.76 0.13 0.19
## ---
## 423: 424 430 1.12 0.02 1.08 1.08
## 424: 425 431 1.12 0.02 1.07 1.08
## 425: 426 432 1.10 0.02 1.06 1.06
## 426: 427 433 1.07 0.02 1.03 1.04
## 427: 428 434 1.08 0.02 1.04 1.05
## Quantile.0.25(R) Median(R) Quantile.0.75(R) Quantile.0.95(R)
## 1: 0.90 1.39 2.04 3.28
## 2: 0.45 0.78 1.25 2.20
## 3: 0.43 0.76 1.21 2.14
## 4: 0.46 0.81 1.30 2.29
## 5: 0.52 0.90 1.44 2.54
## ---
## 423: 1.11 1.12 1.14 1.16
## 424: 1.10 1.12 1.13 1.15
## 425: 1.09 1.10 1.12 1.14
## 426: 1.06 1.07 1.09 1.11
## 427: 1.07 1.08 1.10 1.12
## Quantile.0.975(R) dates Incidence
## 1: 3.76 2020-01-31 0
## 2: 2.58 2020-02-01 0
## 3: 2.51 2020-02-02 0
## 4: 2.69 2020-02-03 0
## 5: 2.99 2020-02-04 0
## ---
## 423: 1.17 2021-03-28 313
## 424: 1.16 2021-03-29 234
## 425: 1.14 2021-03-30 364
## 426: 1.11 2021-03-31 414
## 427: 1.13 2021-04-01 4753. Draw the Curve
Next, create a function that pass the dataset generated by previous step in, and draw it by {echarts4r}.
#' Function for generate a echarts4r object contains a Rt line for one prefecture.
#'
#' @param dataset `data.table` for ploting generated by previous function.
#' @param pref prefecture name.
#'
#' @return a echarts4r object contains a Rt line for one prefecture.
createRtLineFromEstimate <- function(dataset, pref) {
dataset %>%
e_chart(dates, height = "250px") %>%
e_line(
serie = `Mean(R)`,
symbolSize = 0,
name = "Rt"
) %>%
e_bar(
serie = Incidence,
y_index = 1,
barCategoryGap = "10%",
name = "Incidence"
) %>%
e_band(
min = `Quantile.0.05(R)`,
max = `Quantile.0.95(R)`
) %>%
e_tooltip(trigger = "axis") %>%
e_datazoom(
minValueSpan = 28,
startValue = (max(dataset$dates, na.rm = T) - 90)
) %>%
e_legend(top = "10%") %>%
e_y_axis(
name = "Rt",
nameLocation = "middle",
splitLine =
list(lineStyle = list(opacity = 0.4)),
z = 999,
axisLabel = list(inside = T),
axisTick = list(show = F)
) %>%
e_y_axis(
name = "Incidence",
nameGap = 10,
splitLine = list(show = F),
z = 999,
index = 1,
axisTick = list(show = F)
) %>%
e_mark_line(
data = list(yAxis = 1),
symbolSize = 0,
label = list(show = FALSE)
) %>%
e_grid(
left = "8%",
right = "15%"
) %>%
e_title(
text = sprintf(
"%s - Rt: %s±%s（%s）",
pref,
tail(dataset$`Mean(R)`, 1)[[1]],
tail(dataset$`Std(R)`, 1)[[1]],
tail(dataset$dates, 1)[[1]]
),
subtext = sprintf(
"Serial interval: %s±%s",
mean_si, std_si
)
)
}
# Test
createRtLineFromEstimate(tokyo, "東京")# Generate each prefectures Rt line via loop
for (selectedPref in pref) {
dataset <- createEstimateFromSource(byDate, selectedPref)
e <- createRtLineFromEstimate(dataset, selectedPref)
output_list <- c(output_list, list(e))
}
# Render them all at once!
e_arrange(output_list)