实时估计 COVID-19 的 𝑅𝑡（用 R 复现）

请在 Twitter 上关注 Ramnath Vaidyanathan。本分析最初发表在 Google Colab 上。

Kevin Systrom 发布了一篇非常有趣的文章，介绍如何估计 $R_t$，即有效再生数这一指标。在文中，Kevin 将该数值描述为：

在时刻 𝑡，每位传染者平均会感染的人数

Kevin 贴心地将其代码以 notebook 的形式发布在 github 上。我希望在他的模型基础上做更多探索，并将其应用到其他国家。但考虑到我对 R 更为熟练，我首先希望在 R 中复现他的结果。

在阅读本文其余部分前，我强烈建议您先阅读这篇原始文章，以更详细地了解所用的建模方法，因为本文更侧重于代码实现。

建模思路和代码均应归功于 Kevin 👍。本文 notebook 只是我将其模型翻译为 R 的一次拙见尝试。若仍有错误，概由我承担。

Bettencourt & Ribeiro

我们将采用与 Kevin 文中相同的方法，其基础源于 Bettencourt & Ribeiro 的论文。

到达过程建模

第一步是为感染的“到达”过程建模。在统计学中，用于描述到达分布的一个常见选择是泊松分布。因此，若令 $\lambda$ 表示每日平均感染率，则某天出现 $k$ 个新病例的概率为：

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

在此设定下，我们可以为一组 $\lambda$ 构造新病例的概率分布。在开始写代码之前，我们需要加载 tidyverse 包，并自定义一些绘图选项（图形尺寸、绘图风格等）。本节内容对理解主线并非关键，您可跳读，但请务必执行其中的代码。

软件包与选项

# 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))
}
技术支持

── [1mAttaching packages[22m ─────────────────────────────────────── tidyverse 1.3.0 ──

[32m✔[39m [34mggplot2[39m 3.3.0     [32m✔[39m [34mpurrr  [39m 0.3.3
[32m✔[39m [34mtibble [39m 3.0.0     [32m✔[39m [34mdplyr  [39m 0.8.5
[32m✔[39m [34mtidyr  [39m 1.0.2     [32m✔[39m [34mstringr[39m 1.4.0
[32m✔[39m [34mreadr  [39m 1.3.1     [32m✔[39m [34mforcats[39m 0.5.0

── [1mConflicts[22m ────────────────────────────────────────── tidyverse_conflicts() ──
[31m✖[39m [34mdplyr[39m::[32mfilter()[39m masks [34mstats[39m::filter()
[31m✖[39m [34mdplyr[39m::[32mlag()[39m    masks [34mstats[39m::lag()
技术支持

泊松分布

我们可以使用 purrr 中的 crossing 函数，来计算不同 $k$ 与 $\lambda$ 组合下的概率密度。这是一种优雅的方法，但在计算上并非最优，因为它没有充分利用 dpois 的向量化能力。稍后您将看到一种更高效的实现。

# 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)
技术支持

lambda	k	p
10	0	0.0000454
10	1	0.0004540
10	2	0.0022700
10	3	0.0075667
10	4	0.0189166
10	5	0.0378333

我们可以使用 ggplot2 将这些概率可视化。注意如何使用 expression 在图中渲染数学符号。更多信息可参见 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)
  )
技术支持

泊松似然

将到达过程建模为泊松分布，使我们能够用到达率 $\lambda$ 来预测某天新增病例数的分布。但在现实中，我们只观察到到达数量。因此关键问题是：如何从观测到的到达数，反推 $\lambda$ 的分布？幸运的是，答案很简单。

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

$\lambda$ 关于 k 的分布称为似然函数，其数学表达与我们先前使用的概率质量函数相同。我们可以固定观测到的新增病例数（k），并在一系列 $\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
  )
技术支持

连接 $\lambda$ 与 $R_t$

这些都很好，但您可能会问，开头提到的有效感染率 $R_t$ 在哪里？它与到达率 $\lambda$ 有何关系？根据 Bettencourt & Ribeiro 的论文，二者关系相当简单：

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

其中，$\gamma$ 是序列间隔的倒数（COVID-19 约为 4 天），$k_{t-1}$ 是在时间区间 $t-1$ 观察到的新增病例数。

我们可以将该表达式代入 $\lambda$，将似然函数改写为关于 $R_t$ 的形式。

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

在继续之前，我们先设置一些核心参数。

# 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)
技术支持

下面的代码展示了更高效的似然计算方式。我们不再为每种组合分别生成 lambda 与似然，而是利用 R 的向量化与嵌套数据框能力。这使我们能为每一天生成一组 lambda 和似然向量，再借助 tidyr::unnest 将其展开为平面表。

除了计算效率更高外，我认为这段代码更易读，且能以最小的负担贴近数学公式的表达。

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)
技术支持

day	k	r_t	likelihood_r_t
1	40	7.26	0.0000000
1	40	2.54	0.0011292
2	55	3.54	0.0003426
2	55	5.62	0.0000000
1	40	7.03	0.0000000
3	90	10.84	0.0000000

现在我们可以绘制给定观测新增病例数下的条件似然。

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'
  )
技术支持

估计 $R_t$

我们已经很好地将新增病例数 $k$ 与有效感染率 $R_t$ 关联起来。但仍缺少一种可操作的方法，从观测的新增病例时间序列中估计其取值。这时无处不在的贝叶斯法则就派上用场了。Bettencourt & Ribeiro 将 $R_t$ 的概率分布表示为：

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

Kevin 用更朴素的表述解读了该公式：

这意味着，在观察到 𝑘 个新病例后，我们对 𝑅𝑡 的分布的相信等于：

先验 $P(R_t)$，即未使用数据前对 $R_t$ 的信念…

乘以在给定 $R_t$ 下观察到 $k$ 个新病例的似然…

再除以总体上观察到这么多病例的概率。

如果把前一期的后验概率 $P(R_{t-1} | k_{t-1})$ 作为当期的先验 $P(R_t)$，则可将上式改写为：

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

跨期迭代回溯至 $t = 0$，得到：

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

若采用均匀先验 $P(R_0)$，则化简为：

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

我们现在可以用该数学表达式从似然计算后验。同样地，我们用向量化函数 cumprod 一次性生成所有日期的后验，并进行归一化。

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)
技术支持

day	k	r_t	likelihood_r_t	posterior
3	90	5.85	0.0000000	0.0000000
3	90	2.35	0.0033847	0.0025235
3	90	7.68	0.0000000	0.0000000
2	55	0.01	0.0000047	0.0000000
1	40	0.79	0.0000009	0.0000009
3	90	10.81	0.0000000	0.0000000

现在我们可以将后验分布可视化。

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'
  )
技术支持

这些后验现在使我们能够回答一个至关重要的问题：

每天 $R_t$ 的最可能取值是多少？

在估计任何量时，提供估计误差范围非常重要。一个常用方法是使用最高密度区间（HDI）。

另一种总结分布的方法是最高密度区间，即给出覆盖大部分分布（例如 95%）的区间，并保证区间内每一点的可信度都高于区间外任意一点。

在其文章中，Kevin 用了穷举算法来计算 HDI。R 用户更幸运，因为已有现成的 HDInterval 包可用 🎉。

我们现在将估计 $R_t$ 的最可能值及其周围的最高密度区间。由于 HDInterval::hdi 只接受来自分布的随机样本向量，而我们目前只有后验概率，因此需依据后验概率对 $R_t$ 进行随机模拟，以便计算 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)
技术支持

Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
技术支持

day	r_t_most_likely	r_t_lo	r_t_hi
1	3.77	2.45	4.96
2	2.84	2.04	3.65
3	2.90	2.28	3.44

我们可以通过带带状区间的折线图来展示这些估计值及其不确定性。

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
  )
技术支持

尽管 Bettencourt 和 Ribeiro 的方法可行，Kevin 指出了一个实践中的问题（论文中也提到）：由于无法“遗忘”早期高 $R_t$ 的时段，后验可能会卡住并渐近逼近 1。

为解决该问题，Kevin 建议如下做法。

我建议只纳入最近 $m$ 天的似然。这样，算法的先验建立在近期之上，比起整个疫情历程更实用。这个简单但重要的改动带来如下结果：$$ P(R_t|k_t) \propto {\displaystyle \prod^{T}_{t=T-m}}\mathcal{L}\left(k_t|R_t\right) $$

我们将采用该修改，来估计美国各州的 $R_t$。

应用于美国数据

现在是把先前构建的整套方法应用到美国各州 COVID 数据，以估计有效感染率 $R_t$ 的时候了。我们先获取数据，使用纽约时报整理的美国 COVID 数据集。

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

display_head(covid_cases)
技术支持

Parsed with column specification:
cols(
  date = [34mcol_date(format = "")[39m,
  state = [31mcol_character()[39m,
  fips = [31mcol_character()[39m,
  cases = [32mcol_double()[39m,
  deaths = [32mcol_double()[39m
)
技术支持

date	state	fips	cases
2020-01-21	Washington	53	1
2020-01-22	Washington	53	1
2020-01-23	Washington	53	1
2020-01-24	Illinois	17	1
2020-01-24	Washington	53	1
2020-01-25	California	06	1

平滑新增病例

第一步是通过计算每日新增病例并在滚动窗口上进行平滑来准备数据。平滑对于处理报告滞后非常重要，周末的滞后尤为明显。

遵循 Kevin 的方法，我使用 7 天滚动窗口的高斯平滑。

注意我们会将所有数据处理步骤实现为独立函数。虽然这增加了一些开销，但有利于整洁地组合这些步骤，并能方便地在一个或多个州上应用，十分实用。

我们将这些步骤先应用于一个选定州，以便更清晰地查看结果并理解其含义。我选择了纽约州，您也可随意更换为其他州！

# 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()
技术支持

Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
技术支持

state	date	new_cases	new_cases_smooth
New York	2020-03-01	1	1
New York	2020-03-02	0	2
New York	2020-03-03	1	4
New York	2020-03-04	9	9
New York	2020-03-05	11	15
New York	2020-03-06	22	23

我们来绘制 new_cases 的原始值和平滑值，直观展示平滑过程的作用。

🎩 对数据做了变换后及时可视化，始终是良好实践。这有助于建立直觉，也可检验变换是否符合预期！

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()
技术支持

计算似然

第二步是计算似然。我们将沿用之前的方法，但有一处重要不同：我们将计算对数似然而不是似然。这便于在滚动窗口上进行平滑，从而应用 Kevin 提出的修改，即仅使用最近的 $m$ 个区间来计算 $R_t$。

与之前一样，我们会将每天的似然计算为嵌套列表，再用 tidyr::unnest 将结果展开为表格。

🎩 在数据处理中，列表列非常实用。此处使用的设计模式在诸多场景中都十分有用。

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()
技术支持

state	date	new_cases	new_cases_smooth	r_t	likelihood_r_t
New York	2020-03-24	4790	5677	5.36	-4010.9111
New York	2020-03-21	3254	3614	7.33	-5054.4656
New York	2020-04-07	9378	9094	5.90	-10331.2552
New York	2020-03-13	95	118	7.01	-144.3349
New York	2020-04-10	10575	9895	11.69	-107594.6949
New York	2020-04-10	10575	9895	8.46	-35816.2131

计算后验

第三步是计算后验概率。回顾：

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

可改写为对数似然形式：

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

我们可以使用 zoo 包的 rollapplyr 函数，对对数似然进行 7 天滚动求和，再指数化计算后验。最后将后验归一化为 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()
技术支持

state	date	new_cases	new_cases_smooth	r_t	posterior
New York	2020-04-14	7177	8311	0.51	0.0000000
New York	2020-03-19	1770	1921	11.62	0.0000000
New York	2020-03-25	7401	6094	5.46	0.0000000
New York	2020-03-22	4812	4455	0.53	0.0000000
New York	2020-03-20	2950	2753	8.64	0.0000000
New York	2020-03-05	11	15	2.44	0.0020211

我们来可视化所计算的后验概率。注意我们设置了 alpha = 0.2 以减少过度绘制，从而观察到后验的移动。

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()
技术支持

估计 Rt

最后一步是估计 $R_t$ 的取值及其周围的最高密度区间。回顾一下，我们需要根据后验概率对 $r_t$ 进行随机模拟，才能用 HDIInterval::hdi 计算 HDI。

# 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()
技术支持

state	date	r_t_most_likely	r_t_lo	r_t_hi
New York	2020-03-14	2.05	1.68	2.36
New York	2020-03-22	2.33	2.27	2.39
New York	2020-03-07	2.62	1.73	3.46
New York	2020-03-09	2.01	1.32	2.63
New York	2020-04-04	1.18	1.15	1.21
New York	2020-04-06	1.08	1.05	1.11

终于来到“真相时刻”！我们来可视化估计得到的 $R_t$。

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()
技术支持

那么这与 Kevin 文中的估计相比如何？看起来我模型的初始估计与 Kevin 的有些偏差，但这可能源于平滑算法的差异，或初期估计存在较大不确定性。总体曲线形状与峰值位置非常接近。任务完成！🎉

遍历所有州

现在该在所有州上循环并计算这些估计了。我们可以先按 state 分组，将数据拆分为每州一张表，再用 purrr::map_df 为每个州估计 $R_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()
技术支持

estimates_all %>%
  display_random()
技术支持

state	date	r_t_most_likely	r_t_lo	r_t_hi
New Jersey	2020-03-15	2.72	2.01	3.42
Alaska	2020-04-02	0.80	0.04	1.54
Maryland	2020-03-12	1.67	0.00	3.33
Washington	2020-03-10	2.01	1.49	2.52
New Mexico	2020-03-28	1.58	0.89	2.17
South Dakota	2020-04-08	1.63	1.12	2.09

我们可以为所有州的估计绘制小倍数图。有了 ggplot2，这非常简单，只需多加一行代码！

# 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)
技术支持

虽然小倍数图很合适，但基于地理位置查找某个州并不直观。如果能按地理位置把各州排布得更接近呢？

幸运的是，Ryan Hafen 创建了 geofacet 包，让我们仅用一行代码就能把小倍数图布局在地理网格上。

⚠️ 很遗憾，我未能在 colab 上安装 geofacet（我怀疑是 sf 包的系统依赖导致），因此下面的代码块会执行失败。我附上了在本地生成并保存的图像，供您参考其效果。

# ⚠️ 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)
技术支持

最后，我们复刻 Kevin 文中的图：按各州最近一次 $R_t$ 的最可能估计值排序，并根据封锁状态着色。

⚠️ 请注意，这只是描述性图表，无法据此得出有关封锁效果的结论

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)
技术支持

结论与后续工作

本文使用 tidyverse 在 R 中复现了 Kevin 的分析。我们学习了若干计算似然与后验的实用技巧。我在这个过程中收获良多！

我想到不少后续方向，这里重点说一个。鉴于文中的函数仅依赖一个包含三列：location、date、cases 的简单输入表，一个有趣的下一步是将其打包成一个R 包，以便更容易在其他地区估计 $R_t$。

最后，感谢 Kevin Systrom 发布这篇有趣的文章并开源代码。也感谢 DataCamp 的同事们对本文初版提出的精彩反馈。

若您对本文有任何反馈/意见，或想提出扩展建议，欢迎在 Twitter 上联系我 @ramnath_vaidya 或在 github 上找到我。想进一步学习 R 编程

主题

R 数据科学

主题

技术

Estimating COVID-19's 𝑅𝑡 in Real-Time (Replicating in R) Tutorial

使用 AI 探索