Hướng dẫn Ước tính 𝑅𝑡 của COVID-19 theo thời gian thực (Tái lập bằng R)

Hãy theo dõi Ramnath Vaidyanathan trên Twitter. Phân tích này ban đầu được đăng trên Google Colab.

Kevin Systrom đã đăng một bài viết rất thú vị về ước tính $R_t$, thước đo được gọi là số sinh sản hiệu dụng. Trong bài viết, Kevin mô tả con số này là:

Số người bị lây nhiễm bởi mỗi người lây nhiễm tại thời điểm 𝑡

Kevin đã hào phóng công bố mã của mình dưới dạng một notebook trên github. Tôi muốn làm nhiều hơn với mô hình của anh ấy, và cũng áp dụng nó cho các quốc gia khác. Tuy nhiên, vì tôi thành thạo R hơn nhiều, nên trước tiên tôi muốn tái lập kết quả của anh ấy trong R.

Trước khi đọc phần còn lại của bài viết, tôi rất khuyến nghị bạn đọc bài viết gốc để hiểu chi tiết hơn về cách tiếp cận mô hình hoá được sử dụng, vì bài này tập trung nhiều hơn vào phần mã.

Mọi công lao cho cách tiếp cận mô hình và mã đều thuộc về Kevin 👍. Notebook này là nỗ lực nhỏ của tôi nhằm chuyển mô hình của anh ấy sang R. Mọi sai sót còn lại hoàn toàn là của tôi.

Bettencourt & Ribeiro

Chúng ta sẽ dùng cùng cách tiếp cận mà Kevin nêu trong bài viết của anh ấy, dựa trên bài báo của Bettencourt & Ribeiro làm nền tảng.

Mô hình hoá số ca đến

Bước đầu tiên là mô hình hoá quá trình “đến” của các ca nhiễm. Một lựa chọn phổ biến cho phân phối của số ca đến trong giới thống kê là phân phối Poisson. Theo đó, nếu để $\lambda$ biểu thị tốc độ nhiễm trung bình mỗi ngày, thì xác suất chúng ta có thể thấy $k$ ca mới trong một ngày được cho bởi

$$P(k|\lambda) = \frac{\lambda^k e^{-\lambda}}{k!}$$

Với thiết lập này, ta có thể xây dựng phân phối xác suất của số ca mới cho một tập các giá trị $\lambda$. Trước khi viết mã, chúng ta cần nạp bộ gói tidyverse và tuỳ biến một số tuỳ chọn vẽ (kích thước, kiểu biểu đồ, v.v.). Bạn có thể bỏ qua phần này khi đọc, vì nó không trọng yếu với mạch bài, nhưng hãy đảm bảo chạy mã.

Gói & tuỳ chọn

# Load packages
library(tidyverse)

# Plot options

## Jupyter notebooks use the repr package to create viewable representations
## of R objects (https://github.com/IRkernel/repr). I am updating the default
## plot dimensions to 12 x 6.
options(repr.plot.width = 12, repr.plot.height = 6)

## We will use ggplot2 for all plots. I am defining a custom theme here
## that mainly updates the backgrounds and legend position. We set this
## custom theme as the default, and also update the default for line size.
theme_custom <- function(base_size, ...){
  ggplot2::theme_gray(base_size = base_size, ...) +
  ggplot2::theme(
    plot.title = element_text(face = 'bold'),
    plot.subtitle = element_text(color = '#333333'),
    panel.background = element_rect(fill = "#EBF4F7"),
    strip.background = element_rect(fill = "#33AACC"),
    legend.position = "bottom"
  )
}
ggplot2::theme_set(theme_custom(base_size = 20))
ggplot2::update_geom_defaults("line", list(size = 1.5))

# Utility functions

## We will use a utility function to display the head of dataframes.
## Note that we need this hack mainly to add the class 'dataframe' to
## the tables that are printed. This should ideally be handled
## by the `repr` package, and I will be sending a PR.
display_df <- function(x){
  d <- as.character(
    knitr::kable(x, format = 'html', table.attr = "class='dataframe'")
  )
  IRdisplay::display_html(d)
}

display_head <- function(x, n = 6){
   display_df(head(x, n))
}

display_random <- function(x, n = 6){
   display_df(dplyr::sample_n(x, n))
}
Được hỗ trợ bởi

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

✔ ggplot2 3.3.0     ✔ purrr   0.3.3
✔ tibble  3.0.0     ✔ dplyr   0.8.5
✔ tidyr   1.0.2     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats  0.5.0

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
Được hỗ trợ bởi

Phân phối Poisson

Chúng ta có thể dùng hàm crossing từ purrr để tính mật độ xác suất trên mọi kết hợp của $k$ và $\lambda$. Đây là một cách tiếp cận thanh nhã, dù không tối ưu về tính toán, vì ta chưa tận dụng khả năng vector hoá của dpois. Bạn sẽ thấy một cách hiệu quả hơn ở phần sau.

# Number of new cases observed in a day
k = 0:69

# Arrival rate of new infections per day
lambda = c(10, 20, 30, 40)

poisson_densities = crossing(lambda = lambda, k = k) %>%
  mutate(p = dpois(k, lambda))

display_head(poisson_densities)
Được hỗ trợ bởi

Ta có thể trực quan hoá các xác suất này bằng ggplot2. Lưu ý cách chúng ta dùng expression để hiển thị ký hiệu toán học trên biểu đồ. Bạn có thể đọc thêm trong tài liệu của plotmath.

poisson_densities %>%
  # We convert lambda to a factor so that each line gets a discrete color
  mutate(lambda = factor(lambda)) %>%
  ggplot(aes(x = k, y = p, color = lambda)) +
  geom_line() +
  labs(
    title = expression(paste("Probability of k new cases P(k|", lambda, ")")),
    x = 'Number of new cases',
    y = NULL,
    color = expression(lambda)
  )
Được hỗ trợ bởi

Hàm khả năng Poisson

Mô hình hoá quá trình đến bằng phân phối Poisson cho phép ta dự đoán phân phối số ca mới trong một ngày theo hàm của tốc độ đến $\lambda$. Tuy nhiên, trên thực tế, ta chỉ quan sát được số ca đến. Vậy câu hỏi mấu chốt là làm sao đi từ số ca quan sát được để suy ra phân phối của $\lambda$. May mắn là câu trả lời khá đơn giản.

$$L(\lambda| k) = \frac{\lambda^k e^{-\lambda}}{k!}$$

Phân phối của $\lambda$ theo k được gọi là hàm khả năng, và có biểu thức toán học giống với hàm khối xác suất ta dùng trước đó. Ta có thể trực quan hoá hàm khả năng bằng cách cố định số ca mới quan sát được (k), và tính hàm khả năng trên một dải giá trị của $\lambda$.

# Number of new cases observed in a day
k = 20

# Arrival rates of new infections per day
lambdas = seq(1, 45, length = 90)

# Compute likelihood and visualize them
tibble(lambda = lambdas, p = dpois(k, lambdas)) %>%
  ggplot(aes(x = lambda, y = p)) +
  geom_line(color = 'black') +
  labs(
    title = expression(paste("Poisson Likelihood L(", lambda, " | k"[t], ")")),
    x = expression(lambda),
    y = NULL
  )
Được hỗ trợ bởi

Kết nối $\lambda$ và $R_t$

Tất cả những điều trên đều tốt, nhưng bạn có thể đang tự hỏi đâu là hệ số lây nhiễm hiệu dụng, $R_t$, mà ta nhắc đến lúc đầu. Nó liên quan gì đến tốc độ đến $\lambda$. Theo bài báo của Bettencourt & Ribeiro, mối quan hệ khá đơn giản

$$ \lambda = k_{t-1}e^{\gamma(R_t-1)}$$

Lưu ý $\gamma$ ở đây là nghịch đảo của khoảng cách chuỗi lây (khoảng 4 ngày với COVID-19), và $k_{t-1}$ là số ca mới quan sát được trong khoảng thời gian $t-1$.

Ta có thể dùng biểu thức này cho $\lambda$ và viết lại hàm khả năng theo $R_t$.

$$\mathcal{L}\left(R_t|k\right) = \frac{\lambda^k e^{-\lambda}}{k!}$$

Trước khi tiếp tục, hãy thiết lập một số tham số lõi.

# r_t_range is a vector of possible values for R_t
R_T_MAX = 12
r_t_range = seq(0, R_T_MAX, length = R_T_MAX*100 + 1)

# Gamma is 1/serial interval
# https://wwwnc.cdc.gov/eid/article/26/6/20-0357_article
GAMMA = 1/4

# New cases by day
k =  c(20, 40, 55, 90)
Được hỗ trợ bởi

Đoạn mã dưới đây cho thấy cách tính hàm khả năng hiệu quả hơn. Thay vì tạo lambdas và khả năng cho từng kết hợp riêng lẻ, ta tận dụng tính vector hoá của các hàm trong R và khả năng có dataframe lồng nhau. Cách này cho phép tạo một vector lambdas và likelihoods cho mỗi ngày, rồi dùng tidyr::unnest để bung thành bảng phẳng.

Ngoài hiệu quả tính toán, tôi cho rằng mã này cũng dễ đọc hơn và chuyển tải công thức toán học với chi phí tối thiểu.

likelihoods <- tibble(day = seq_along(k) - 1, k = k) %>%
  # Compute a vector of likelihoods
  mutate(
    r_t = list(r_t_range),
    lambda = map(lag(k, 1), ~ .x * exp(GAMMA * (r_t_range - 1))),
    likelihood_r_t = map2(k, lambda, ~ dpois(.x, .y)/sum(dpois(.x, .y)))
  ) %>%
  # Ignore the 0th day
  filter(day > 0) %>%
  # Unnest the data to flatten it.
  select(-lambda) %>%
  unnest(c(r_t, likelihood_r_t))

display_random(likelihoods)
Được hỗ trợ bởi

Bây giờ chúng ta có thể vẽ hàm khả năng có điều kiện theo số ca mới quan sát được.

likelihoods %>%  
  ggplot(aes(x = r_t, y = likelihood_r_t, color = factor(k))) +
  geom_line() +
  labs(
    title = expression(paste("Likelihood of R"[t], " given k")),
    subtitle = expression(paste("L(R"[t], "|k)")),
    x = expression("R"[t]),
    y = NULL, color = 'k'
  )
Được hỗ trợ bởi

Ước tính $R_t$

Chúng ta đã tiến khá xa trong việc liên hệ số ca mới $k$ và hệ số lây nhiễm hiệu dụng $R_t$. Tuy vậy, ta vẫn chưa có cách cụ thể để ước tính các giá trị của nó từ chuỗi thời gian số ca mới quan sát được. Đây là lúc định lý Bayes phát huy tác dụng. Bettencourt & Ribeiro biểu diễn phân phối xác suất của $R_t$ như sau

$$ P(R_t|k_t)=\frac{P(R_t)\cdot\mathcal{L}(k_t|R_t)}{P(k_t)} $$

Kevin diễn giải phương trình phức tạp này theo cách đơn giản hơn nhiều.

Điều này nói rằng, sau khi thấy 𝑘 ca mới, chúng ta cho rằng phân phối của 𝑅𝑡 bằng:

• tiên nghiệm P(R_t) về giá trị của nó khi chưa có dữ liệu ...
• nhân với khả năng thấy $k$ ca mới khi biết $R_t$ ...
• chia cho xác suất nói chung của việc thấy chừng đó ca.

Nếu ta dùng xác suất hậu nghiệm của kỳ trước, $P(R_{t-1} | k_{t-1})$, làm tiên nghiệm, $P(R_t)$, cho kỳ hiện tại, thì ta có thể viết lại phương trình trước thành:

$$ P(R_t|k_t) \propto P(R_{t-1}|k_{t-1})\cdot\mathcal{L}(k_t|R_t) $$

Lặp qua các kỳ ngược về $t = 0$, ta được

$$ P(R_t|k_t) \propto P(R_0) \cdot {\displaystyle \prod^{T}_{t=0}}\mathcal{L}(k_t|R_t) $$

Với tiên nghiệm đều $P(R_0)$, điều này rút gọn thành:

$$ P(R_t|k_t) \propto {\displaystyle \prod^{T}_{t=0}}\mathcal{L}\left(k_t|R_t\right) $$

Giờ đây chúng ta có thể dùng biểu thức này để tính hậu nghiệm từ các khả năng. Một lần nữa, ta dùng hàm vector hoá cumprod để tạo hậu nghiệm cho tất cả các ngày cùng lúc, rồi chuẩn hoá.

posteriors <- likelihoods %>%
  group_by(r_t) %>%
  arrange(day) %>%
  mutate(posterior = cumprod(likelihood_r_t)) %>%
  group_by(k) %>%
  mutate(posterior = posterior / sum(posterior)) %>%
  ungroup()

display_random(posteriors)
Được hỗ trợ bởi

Giờ chúng ta có thể trực quan hoá phân phối hậu nghiệm.

posteriors %>%
  ggplot(aes(x = r_t, y = posterior, color = factor(day))) +
  geom_line() +
  labs(
    title = expression(paste("Posterior probability of R"[t], " given k")),
    subtitle = expression(paste("P(R"[t], "| k)")),
    x = expression("R"[t]), y = NULL, color = 'day'
  )
Được hỗ trợ bởi

Các hậu nghiệm này cho phép chúng ta trả lời câu hỏi quan trọng nhất:

Giá trị khả dĩ nhất của $R_t$ mỗi ngày là bao nhiêu?

Khi ước tính bất kỳ đại lượng nào, điều rất quan trọng là đưa ra cảm nhận về sai số xung quanh ước tính đó. Một cách phổ biến là sử dụng khoảng mật độ cao nhất (HDI).

Một cách khác để tóm tắt một phân phối là khoảng mật độ cao nhất, tóm tắt bằng một khoảng bao quát phần lớn phân phối, chẳng hạn 95% của nó, sao cho mọi điểm bên trong khoảng có độ tin cậy cao hơn bất kỳ điểm nào bên ngoài khoảng.

Trong bài viết, Kevin triển khai một thuật toán vét cạn để tính HDI. Người dùng R may mắn hơn vì đã có sẵn gói HDInterval với hiện thực 🎉.

Chúng ta sẽ ước tính giá trị khả dĩ nhất của $R_t$, và khoảng mật độ cao nhất xung quanh nó. Vì hàm HDInterval::hdi chỉ hoạt động với một vector các giá trị ngẫu nhiên từ phân phối, trong khi ta mới ước tính các xác suất hậu nghiệm, nên ta cần mô phỏng các giá trị ngẫu nhiên của $R_t$ theo các xác suất hậu nghiệm để tính HDI.

# Install and load HDInterval package
install.packages("HDInterval")
library(HDInterval)

# Compute the most likely value of r_t and the highest-density interval
estimates <- posteriors %>%
  group_by(day) %>%
  summarize(
    r_t_simulated = list(sample(r_t_range, 10000, replace = TRUE, prob = posterior)),
    r_t_most_likely = r_t_range[which.max(posterior)]
  ) %>%
  mutate(
    r_t_lo = map_dbl(r_t_simulated, ~ hdi(.x)[1]),
    r_t_hi = map_dbl(r_t_simulated, ~ hdi(.x)[2])
  ) %>%
  select(-r_t_simulated)

display_head(estimates)
Được hỗ trợ bởi

Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
Được hỗ trợ bởi

Ta có thể trực quan hoá các ước tính này và mức độ bất định xung quanh chúng bằng biểu đồ đường kèm dải băng.

estimates %>%
  ggplot(aes(x = day, y = r_t_most_likely)) +
  geom_point(color = "#ffc844", size = 5) +
  geom_line(color = 'black') +
  geom_ribbon(aes(ymin = r_t_lo, ymax = r_t_hi), fill = "#ffc844", alpha = 0.3) +
  labs(
    title = expression(paste('R'[t], ' by day')),
    subtitle = "The band represents the highest density interval",
    x = 'Day', y = NULL
  )
Được hỗ trợ bởi

Mặc dù cách tiếp cận của Bettencourt và Ribeiro hoạt động ổn, Kevin chỉ ra một vấn đề thực tế (cũng được nêu trong bài báo), rằng hậu nghiệm có thể bị kẹt và tiệm cận về 1, do không thể “quên” những thời kỳ có $R_t$ cao.

Để khắc phục, Kevin gợi ý cách làm sau.

Tôi đề xuất chỉ gộp $m$ ngày cuối của hàm khả năng. Bằng cách này, tiên nghiệm của thuật toán được xây dựng dựa trên quá khứ gần, hữu ích hơn nhiều so với toàn bộ lịch sử dịch. Thay đổi đơn giản nhưng quan trọng này dẫn đến: $$ P(R_t|k_t) \propto {\displaystyle \prod^{T}_{t=T-m}}\mathcal{L}\left(k_t|R_t\right) $$

Chúng ta sẽ dùng điều chỉnh này để ước tính $R_t$ trên các bang ở Mỹ.

Ứng dụng với dữ liệu Mỹ

Đến lúc mang toàn bộ bộ công cụ này áp dụng lên dữ liệu COVID từ các bang của Mỹ để ước tính hệ số lây nhiễm hiệu dụng $R_t$. Bắt đầu bằng việc lấy dữ liệu. Chúng ta sẽ dùng bộ dữ liệu COVID cho Mỹ do New York Times biên soạn.

url = 'https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv'
covid_cases <- readr::read_csv(url)

display_head(covid_cases)
Được hỗ trợ bởi

Parsed with column specification:
cols(
  date = col_date(format = ""),
  state = col_character(),
  fips = col_character(),
  cases = col_double(),
  deaths = col_double()
)
Được hỗ trợ bởi

Làm trơn số ca mới

Bước đầu là chuẩn bị dữ liệu bằng cách tính số ca mới mỗi ngày và làm trơn chúng trên một cửa sổ trượt. Làm trơn là cần thiết để tính đến độ trễ báo cáo. Độ trễ này đặc biệt rõ vào cuối tuần.

Theo cách tiếp cận của Kevin, tôi dùng bộ làm trơn Gaussian với cửa sổ 7 ngày.

Lưu ý chúng ta sẽ hiện thực mọi bước xử lý dữ liệu dưới dạng hàm độc lập. Dù có thêm chút chi phí, điều này cho phép ghép các bước gọn gàng, và áp dụng cho một hay nhiều bang — rất hữu ích.

Chúng ta sẽ áp dụng các bước này cho một bang được chọn, để dễ xem kết quả và hiểu chúng. Tôi chọn New York, nhưng bạn có thể đổi sang bang mình muốn!

# Install the smoother package
install.packages("smoother")

#' Compute new cases and smooth them
smooth_new_cases <- function(cases){
  cases %>%
    arrange(date) %>%
    mutate(new_cases = c(cases[1], diff(cases))) %>%
    mutate(new_cases_smooth = round(
      smoother::smth(new_cases, window = 7, tails = TRUE)
    )) %>%
    select(state, date, new_cases, new_cases_smooth)
}

state_selected <- "New York"
covid_cases %>%
  filter(state == state_selected) %>%
  smooth_new_cases() %>%
  display_head()
Được hỗ trợ bởi

Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
Được hỗ trợ bởi

Hãy vẽ cả giá trị thô lẫn giá trị đã làm trơn của new_cases để hình dung quá trình làm trơn đang làm gì.

🎩 Luôn là thực hành tốt khi trực quan hoá kết quả sau mỗi bước biến đổi dữ liệu. Điều này giúp bạn xây dựng trực giác và kiểm tra nhanh xem biến đổi có đúng như dự định không!

plot_new_cases <- function(cases){
  cases %>%
    ggplot(aes(x = date, y = new_cases)) +
    geom_line(linetype = 'dotted', color = 'gray40') +
    geom_line(aes(y = new_cases_smooth), color = "#14243e") +
    labs(
      title = "New cases per day",
      subtitle = unique(cases$state),
      x = NULL, y = NULL
    )
}

covid_cases %>%
  filter(state == state_selected) %>%
  smooth_new_cases() %>%
  plot_new_cases()
Được hỗ trợ bởi

Tính hàm khả năng

Bước thứ hai là tính các likelihood. Ta sẽ dùng lại cách trước đó, nhưng có một khác biệt đáng chú ý. Ta sẽ tính log-likelihood thay vì likelihood. Điều này giúp dễ làm trơn trên cửa sổ trượt để áp dụng điều chỉnh của Kevin: chỉ dùng $m$ khoảng gần nhất để tính $R_t$.

Như trước, ta tính likelihood cho từng ngày dưới dạng danh sách lồng và sau đó bung kết quả thành bảng bằng tidyr::unnest.

🎩 Các cột dạng danh sách cực kỳ hữu ích khi làm việc với dữ liệu. Mẫu thiết kế dùng ở đây rất hữu dụng trong nhiều tình huống.

compute_likelihood <- function(cases){
  likelihood <- cases %>%
    filter(new_cases_smooth > 0) %>%
    mutate(
      r_t = list(r_t_range),
      lambda = map(lag(new_cases_smooth, 1), ~ .x * exp(GAMMA * (r_t_range - 1))),
      likelihood_r_t = map2(new_cases_smooth, lambda, dpois, log = TRUE)
    ) %>%
    slice(-1) %>%
    select(-lambda) %>%
    unnest(c(likelihood_r_t, r_t))
}

covid_cases %>%
  filter(state == state_selected) %>%
  smooth_new_cases() %>%
  compute_likelihood() %>%
  display_random()
Được hỗ trợ bởi

Tính hậu nghiệm

Bước ba là tính các xác suất hậu nghiệm. Nhớ rằng

$$ P(R_t|k_t) \propto {\displaystyle \prod^{T}_{t=T-m}}\mathcal{L}\left(k_t|R_t\right) $$

Ta có thể viết lại theo log-likelihood

$$ P(R_t|k_t) \propto {\displaystyle \exp \big( \sum^{T}_{t=T-m}}\log(\mathcal{L}\left(k_t|R_t\right)) \big) $$

Ta dùng hàm rollapplyr từ gói zoo để tính tổng cuộn 7 ngày của log-likelihood, rồi lấy mũ để tính hậu nghiệm. Cuối cùng, chuẩn hoá hậu nghiệm về 1.0.

compute_posterior <- function(likelihood){
  likelihood %>%
    arrange(date) %>%
    group_by(r_t) %>%
    mutate(posterior = exp(
      zoo::rollapplyr(likelihood_r_t, 7, sum, partial = TRUE)
    )) %>%
    group_by(date) %>%
    mutate(posterior = posterior / sum(posterior, na.rm = TRUE)) %>%
    # HACK: NaNs in the posterior create issues later on. So we remove them.
    mutate(posterior = ifelse(is.nan(posterior), 0, posterior)) %>%
    ungroup() %>%
    select(-likelihood_r_t)
}

covid_cases %>%
  filter(state == state_selected) %>%
  smooth_new_cases() %>%
  compute_likelihood() %>%
  compute_posterior() %>%
  display_random()
Được hỗ trợ bởi

Hãy trực quan hoá các xác suất hậu nghiệm vừa tính. Lưu ý cách ta đặt alpha = 0.2 để giảm chồng lấn nét vẽ. Điều này giúp thấy rõ sự dịch chuyển của hậu nghiệm.

plot_posteriors <- function(posteriors){
  posteriors %>%
    ggplot(aes(x = r_t, y = posterior, group = date)) +
    geom_line(alpha = 0.2) +
    labs(
      title = expression(paste("Daily Posterior of R"[t], " by day")),
      subtitle = unique(posteriors$state),
      x = '',
      y = ''
    ) +
    coord_cartesian(xlim = c(0.4, 4)) +
    theme(legend.position = 'none')
}

covid_cases %>%
  filter(state == state_selected) %>%
  smooth_new_cases() %>%
  compute_likelihood() %>%
  compute_posterior() %>%
  plot_posteriors()
Được hỗ trợ bởi

Ước tính Rt

Bước cuối cùng là ước tính giá trị của $R_t$ và khoảng mật độ cao nhất xung quanh chúng. Nhớ rằng ta cần mô phỏng các giá trị ngẫu nhiên cho $r_t$ theo xác suất hậu nghiệm để áp dụng hàm HDIInterval::hdi tính khoảng mật độ cao nhất.

# Estimate R_t and a 95% highest-density interval around it
estimate_rt <- function(posteriors){
  posteriors %>%
    group_by(state, date) %>%
    summarize(
      r_t_simulated = list(sample(r_t_range, 10000, replace = TRUE, prob = posterior)),
      r_t_most_likely = r_t_range[which.max(posterior)]
    ) %>%
    mutate(
      r_t_lo = map_dbl(r_t_simulated, ~ hdi(.x)[1]),
      r_t_hi = map_dbl(r_t_simulated, ~ hdi(.x)[2])
    ) %>%
    select(-r_t_simulated)
}

covid_cases %>%
  filter(state == state_selected) %>%
  smooth_new_cases() %>%
  compute_likelihood() %>%
  compute_posterior() %>%
  estimate_rt() %>%
  display_random()
Được hỗ trợ bởi

Cuối cùng cũng đến giây phút quyết định! Hãy trực quan hoá các giá trị $R_t$ ước tính

plot_estimates <- function(estimates){
  estimates %>%
    ggplot(aes(x = date, y = r_t_most_likely)) +
    geom_point(color = "darkorange", alpha = 0.8, size = 4) +
    geom_line(color = "#14243e") +
    geom_hline(yintercept = 1, linetype = 'dashed') +
    geom_ribbon(
      aes(ymin = r_t_lo, ymax = r_t_hi),
      fill = 'darkred',
      alpha = 0.2
    ) +
    labs(
      title = expression('Real time R'[t]), x = '', y = '',
      subtitle = unique(estimates$state)
    ) +
    coord_cartesian(ylim = c(0, 4))
}

covid_cases %>%
  filter(state == state_selected) %>%
  smooth_new_cases() %>%
  compute_likelihood() %>%
  compute_posterior() %>%
  estimate_rt() %>%
  plot_estimates()
Được hỗ trợ bởi

Vậy điều này so với các ước tính trong bài của Kevin thế nào? Có vẻ như các ước tính ban đầu trong mô hình của tôi lệch so với của Kevin. Tuy nhiên, đó có thể do khác biệt trong thuật toán làm trơn, hoặc do bất định lớn quanh các ước tính ban đầu. Hình dáng tổng thể của đường cong và các đỉnh trông khá khớp. Hoàn thành nhiệm vụ! 🎉

Lặp qua tất cả các bang

Giờ là lúc lặp qua tất cả các bang và tính các ước tính này. Ta có thể dễ dàng làm điều đó bằng cách nhóm theo state, tách dữ liệu thành mỗi bang một bảng, rồi dùng purrr::map_df để ước tính $R_t$ cho từng bang và gộp lại thành một bảng duy nhất.

# ⚠️This function can take a couple of minutes to run
#   as it loops across all states
estimates_all <- covid_cases %>%
  filter(date >= "2020-03-01") %>%
  group_by(state) %>%
  # Ignore states that have not reached 100 infections
  filter(max(cases) > 100 ) %>%
  group_split() %>%
  map_df(~ {
    .x %>%
      smooth_new_cases() %>%
      compute_likelihood() %>%
      compute_posterior() %>%
      estimate_rt()
  }) %>%
  ungroup()
Được hỗ trợ bởi

estimates_all %>%
  display_random()
Được hỗ trợ bởi

Giờ ta có thể tạo biểu đồ small multiples các ước tính cho tất cả các bang. Việc dùng ggplot2 khiến điều này trở nên rất dễ — chúng ta chỉ cần thêm một dòng mã!

# Increase plot height and width
options(repr.plot.height = 40, repr.plot.width = 20)
estimates_all %>%
  plot_estimates() +
  facet_wrap(~ state, ncol = 4) +
  labs(subtitle = "")

# Reset plot dimensions
options(repr.plot.height = 12, repr.plot.width = 8)
Được hỗ trợ bởi

Mặc dù biểu đồ small multiples rất ổn, nhưng khó để tìm một bang theo vị trí địa lý. Nếu ta có thể vẽ mỗi bang gần vị trí địa lý của nó thì sao?

May thay, Ryan Hafen đã tạo gói geofacet cho phép sắp xếp các small multiples trên một lưới địa lý chỉ với một dòng mã.

⚠️ Tiếc là tôi không cài được gói geofacet trên colab (tôi nghi do phụ thuộc hệ thống của gói sf), nên khối mã dưới sẽ chạy lỗi. Tôi đã đính kèm bản lưu biểu đồ tạo cục bộ để bạn hình dung.

# ⚠️ CAUTION: This code will error on colab
# options(repr.plot.height = 40, repr.plot.width = 20)
# estimates_all %>%
#  mutate(state = state.abb[match(state, state.name)]) %>%
#  plot_estimates() +
#  geofacet::facet_geo(~ state, ncol = 4) +
#  labs(subtitle = "") +
#  theme(strip.text = element_text(hjust = 0))
# options(repr.plot.height = 12, repr.plot.width = 8)
Được hỗ trợ bởi

Cuối cùng, hãy tái tạo biểu đồ trong bài của Kevin, sắp xếp các bang theo giá trị $R_t$ khả dĩ nhất gần đây, và tô màu theo trạng thái phong toả.

⚠️ Xin lưu ý đây chỉ là biểu đồ mô tả và không thể rút ra kết luận về hiệu quả của phong toả

options(repr.plot.width = 20, repr.plot.height = 8)
no_lockdown = c('North Dakota', 'South Dakota', 'Nebraska', 'Iowa', 'Arkansas')
partial_lockdown = c('Utah', 'Wyoming', 'Oklahoma')
estimates_all %>%
  group_by(state) %>%
  filter(date == max(date)) %>%
  ungroup() %>%
  mutate(state = forcats::fct_reorder(state, r_t_most_likely)) %>%
  mutate(lockdown = case_when(
    state %in% no_lockdown ~ 'None',
    state %in% partial_lockdown ~ 'Partial',
    TRUE ~ "Full"
  )) %>%
  ggplot(aes(x = state, y = r_t_most_likely)) +
  geom_col(aes(fill = lockdown)) +
  geom_hline(yintercept = 1, linetype = 'dotted') +
  geom_errorbar(aes(ymin = r_t_lo, ymax = r_t_hi), width = 0.2) +
  scale_fill_manual(values = c(None = 'darkred', Partial = 'gray50', Full = 'gray70')) +
  labs(
    title = expression(paste("Most Recent R"[t], " by state")),
    x = '', y = ''
  ) +
  theme(axis.text.x.bottom = element_text(angle = 90, hjust = 1, vjust = 0.5))
options(repr.plot.width = 12, repr.plot.height = 5)
Được hỗ trợ bởi

Kết luận và Bước tiếp theo

Trong bài viết này, chúng ta đã tái lập phân tích của Kevin trong R, dùng tidyverse. Chúng ta học được vài mẹo hữu ích để tính likelihood và hậu nghiệm. Tôi học được rất nhiều trong quá trình này!

Tôi có thể nghĩ ra vài bước tiếp theo, nhưng xin nhấn mạnh một ý lớn. Vì các hàm trong bài dựa vào một bảng đầu vào khá đơn giản với ba cột: location, date, cases, một bước thú vị tiếp theo là biến nó thành một gói R. Điều này sẽ giúp ước tính $R_t$ cho các khu vực khác dễ dàng hơn.

Tôi muốn kết thúc bằng lời cảm ơn Kevin Systrom vì đã đăng một bài viết thú vị và chia sẻ mã nguồn. Tôi cũng cảm ơn các đồng nghiệp tại DataCamp vì phản hồi xuất sắc cho phiên bản đầu tiên của bài này.

Nếu bạn có phản hồi/bình luận về bài viết, hoặc gợi ý mở rộng khác, bạn có thể tìm tôi trên twitter @ramnath_vaidya và github.  Để tìm hiểu thêm về lập trình R