COVID-19'un 𝑅𝑡 Değerini Gerçek Zamanlı Tahmin Etme (R ile Yeniden Üretim) Eğitimi

Twitter'da Ramnath Vaidyanathan'ı takip edin. Bu analiz ilk olarak Google Colab üzerinde yayımlandı.

Kevin Systrom, $R_t$ (etkili üreme sayısı) tahmini üzerine gerçekten ilginç bir yazı yayımladı. Kevin, yazısında bu sayıyı şöyle tanımlıyor:

t anında, bulaştırıcı bir kişi başına enfekte olan kişi sayısı

Kevin, kodunu GitHub'da bir notebook olarak paylaşma nezaketinde bulundu. Modeliyle daha fazlasını yapmak ve onu diğer ülkelere de uygulamak istedim. Ancak, R konusunda çok daha yetkin olduğum için, önce sonuçlarını R'de yeniden üretmek istedim.

Bu yazının kalanını okumadan önce, kullanılan modelleme yaklaşımını daha ayrıntılı anlamak için orijinal makaleyi okumanızı kuvvetle öneririm; çünkü bu yazı daha çok koda odaklanıyor.

Modelleme yaklaşımı ve kod için tüm kredi Kevin'a ait 👍. Bu notebook, onun modelini R'ye çevirmeye mütevazı bir girişimdir. Kalan tüm hatalar tamamen bana aittir.

Bettencourt & Ribeiro

Kevin'ın yazısında özetlediği ve temelini Bettencourt & Ribeiro makalesinden alan aynı yaklaşımı kullanacağız.

Gelişleri modelleme

İlk adım, enfeksiyonların 'geliş' sürecini modellemektir. İstatistikçiler arasında gelişlerin dağılımı için popüler bir seçim Poisson Dağılımıdır. Buna göre, eğer $\lambda$ günlük ortalama enfeksiyon hızını temsil ederse, bir günde $k$ yeni vaka görme olasılığı şu şekilde verilir

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

Bu kurulumla, bir $\lambda$ kümesi için yeni vakaların olasılık dağılımını oluşturabiliriz. Koda başlamadan önce, tidyverse paket setini yüklememiz ve bazı çizim seçeneklerini (grafik boyutları, çizim vb.) özelleştirmemiz gerekiyor. Bu bölümü okumayı atlayabilirsiniz; anlatım açısından kritik değil, ancak kodu çalıştırdığınızdan emin olun.

Paketler & seçenekler

# 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))
}
Altyapı:

── 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()
Altyapı:

Poisson dağılımı

purrr içindeki crossing fonksiyonunu kullanarak, $k$ ve $\lambda$ kombinasyonları boyunca olasılık yoğunluklarını hesaplayabiliriz. Bu, dpois'in vektörleştirme yeteneklerinden tam olarak yararlanmadığımız için hesaplama açısından optimal olmasa da şık bir yaklaşımdır. Daha verimli bir yaklaşımı ileride göreceksiniz.

# 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)
Altyapı:

Bu olasılıkları ggplot2 ile görselleştirebiliriz. Grafiklerde matematiksel sembolleri göstermek için expression kullanımına dikkat edin. Bununla ilgili daha fazla bilgiyi plotmath dokümantasyonunda bulabilirsiniz.

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)
  )
Altyapı:

Poisson olabilirlik

Geliş sürecini Poisson dağılımı olarak modellemek, bir günde yeni vaka dağılımını geliş hızı $\lambda$'nın bir fonksiyonu olarak tahmin etmemizi sağlar. Ancak gerçekte yalnızca geliş sayılarını gözlemleriz. Dolayısıyla asıl soru, gözlemlenen geliş sayısından $\lambda$ dağılımını nasıl çıkarabileceğimizdir. Neyse ki yanıt basit.

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

$\lambda$'nın k üzerindeki dağılımına olabilirlik fonksiyonu denir ve önce kullandığımız olasılık kütle fonksiyonu ile aynı matematiksel ifadeye sahiptir. Gözlemlenen yeni vaka sayısını (k) sabitleyip $\lambda$ için bir değer aralığında olabilirlik fonksiyonunu hesaplayarak görselleştirebiliriz.

# 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
  )
Altyapı:

$\lambda$ ile $R_t$ arasındaki bağlantı

Tüm bunlar güzel, ancak başta bahsettiğimiz etkili bulaşma hızı $R_t$ nerede diye merak ediyor olabilirsiniz. $\lambda$ geliş hızıyla nasıl ilişkilidir? Bettencourt & Ribeiro'nun makalesine göre ilişki oldukça basit

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

Burada $\gamma$'nın seri aralığının tersi olduğuna dikkat edin (COVID-19 için yaklaşık 4 gün) ve $k_{t-1}$, $t-1$ zaman aralığında gözlemlenen yeni vaka sayısıdır.

Bu $\lambda$ ifadesini kullanarak olabilirlik fonksiyonunu $R_t$ cinsinden yeniden yazabiliriz.

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

Devam etmeden önce bazı temel parametreleri ayarlayalım.

# 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)
Altyapı:

Aşağıdaki kod, olabilirlikleri hesaplamanın daha verimli bir yolunu gösteriyor. Her kombinasyon için lambda ve olabilirlikleri ayrı ayrı üretmek yerine, R fonksiyonlarının vektörleştirilmiş doğasını ve iç içe veri çerçevelerine sahip olma yeteneğini kullanıyoruz. Bu, her gün için birer vektör halinde lambda ve olabilirlik üretmemizi ve ardından tidyr::unnest ile düz bir tabloya açmamızı sağlıyor.

Hesaplama verimliliğinin ötesinde, bu kodun daha okunaklı olduğunu ve matematiksel formülleri asgari ek yükle çevirdiğini söyleyebilirim.

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)
Altyapı:

Artık gözlemlenen yeni vaka sayısına koşullu olabilirliği çizebiliriz.

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'
  )
Altyapı:

$R_t$ tahmini

Yeni vaka sayısı $k$ ile etkili bulaşma hızı $R_t$ arasındaki bağlantıda iyi ilerleme kaydettik. Ancak hâlâ, gözlemlenen yeni vaka zaman serisinden değerleri tahmin etmeye yönelik somut bir yaklaşımımız yok. İşte burada her yerde karşımıza çıkan Bayes kuralı son derece kullanışlıdır. Bettencourt & Ribeiro, $R_t$'nin olasılık dağılımını şöyle ifade eder:

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

Kevin bu karmaşık denklemi çok daha basit terimlerle açıyor.

Bu, 𝑘 yeni vaka gördüğümüzde, 𝑅𝑡'nin dağılımına ilişkin inancımızın şu eşitliğe sahip olduğunu söyler:

• Veri olmadan $P(R_t)$'nin değerine ilişkin öncül inançlar ...
• $R_t$ verildiğinde $k$ yeni vaka görmenin olabilirliği ile çarpılır ...
• genel olarak bu kadar çok vaka görme olasılığına bölünür.

Önceki dönemin artgönderimini, yani $P(R_{t-1} | k_{t-1})$'i, cari dönemin öncülü $P(R_t)$ olarak kullanırsak, önceki denklemi şöyle yeniden yazabiliriz:

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

Dönemler arasında $t = 0$'a kadar yinelediğimizde elde ederiz

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

Uniform bir öncül $P(R_0)$ ile bu şuna indirgenir:

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

Bu matematiksel ifadeyi kullanarak olabilirliklerden artgönderimi hesaplayabiliriz. Yine, tüm günler için artgönderimleri bir kerede üretmek üzere vektörleştirilmiş cumprod fonksiyonunu kullanır ve ardından normalleştiririz.

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)
Altyapı:

Artık artgönderim dağılımını görselleştirebiliriz.

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'
  )
Altyapı:

Bu artgönderimler şimdi en kritik soruyu yanıtlamamızı sağlıyor:

Her gün için $R_t$'nin en olası değeri nedir?

Herhangi bir niceliği tahmin ederken, tahmini çevreleyen hataya dair bir fikir vermek çok önemlidir. Bunun popüler bir yolu en yüksek yoğunluk aralıklarını kullanmaktır.

Bir dağılımı özetlemenin bir başka yolu da, dağılımın örneğin %95'ini kapsayan ve aralık içindeki her noktanın, aralık dışındaki herhangi bir noktadan daha yüksek güvenilirliğe sahip olduğu bir aralığı belirten en yüksek yoğunluk aralığıdır.

Kevin yazısında HDI'yi hesaplamak için kaba kuvvet bir algoritma uygulamış. R kullanıcıları daha şanslı; bunu gerçekleştiren HDInterval adlı bir paket zaten var 🎉.

Şimdi $R_t$'nin en olası değerini ve onu çevreleyen en yüksek yoğunluk aralığını tahmin edeceğiz. HDInterval::hdi fonksiyonu yalnızca bir dağılımdan rastgele değerler vektörü ile çalıştığından ve biz yalnızca artgönderim olasılıklarını tahmin ettiğimizden, en yüksek yoğunluk aralıklarını hesaplamak için artgönderim olasılıklarını kullanarak $R_t$'nin rastgele değerlerini simüle etmemiz gerekecek.

# 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)
Altyapı:

Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
Altyapı:

Bu tahminleri ve çevrelerindeki belirsizliği şeritli bir çizgi grafikle görselleştirebiliriz.

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
  )
Altyapı:

Bettencourt ve Ribeiro'nun yaklaşımı işe yarasa da Kevin pratik bir soruna dikkat çekiyor (makalede de belirtilmiş): artgönderimin takılı kalıp asimptotik olarak 1'e yaklaşabilmesi. Bunun nedeni, yüksek $R_t$ dönemlerini unutamamasıdır.

Bu sorunu gidermek için Kevin aşağıdaki yaklaşımı öneriyor.

Bu nedenle, olabilirlik fonksiyonunun yalnızca son $m$ gününü dikkate almayı öneriyorum. Bunu yaparak, algoritmanın öncülü yakın geçmişe dayanır; bu da salgının tüm geçmişine göre çok daha faydalı bir öncül olur. Bu basit ama önemli değişiklik şu sonucu verir: $$ P(R_t|k_t) \propto {\displaystyle \prod^{T}_{t=T-m}}\mathcal{L}\left(k_t|R_t\right) $$

ABD'deki farklı eyaletlerde $R_t$ değerlerini tahmin etmek için bu değişikliği kullanacağız.

ABD verilerine uygulama

Artık kurduğumuz bu yapıyı alıp ABD eyaletlerinden COVID verilerine uygulayarak etkili bulaşma hızı $R_t$'yi tahmin etme zamanı. Verileri çekerek başlayalım. New York Times tarafından oluşturulan ABD'ye ait COVID veri setini kullanacağız.

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

display_head(covid_cases)
Altyapı:

Parsed with column specification:
cols(
  date = col_date(format = ""),
  state = col_character(),
  fips = col_character(),
  cases = col_double(),
  deaths = col_double()
)
Altyapı:

Yeni vakaları düzleştirme

İlk adım, her günkü yeni vaka sayısını hesaplayıp bir hareketli pencerede düzleştirerek veriyi hazırlamaktır. Düzleştirme, raporlama gecikmelerini hesaba katmak için gereklidir. Bu gecikmenin özellikle hafta sonlarında belirgin olduğu görülmüştür.

Kevin'ın yaklaşımını izleyerek, 7 günlük hareketli pencereye sahip bir Gaussian düzleştirici kullanıyorum.

Tüm veri işleme adımlarını bağımsız fonksiyonlar olarak uygulayacağımızı unutmayın. Bu biraz ek yük getirse de adımları temizce birleştirmemizi ve bir veya daha fazla eyalete uygulamamızı sağlar; bu da oldukça işe yarayacaktır.

Sonuçları görmek ve anlamlandırmak daha kolay olsun diye bu adımları seçili bir eyalete uygulayacağız. Ben New York'u seçtim, ancak siz istediğiniz bir eyaletle değiştirebilirsiniz!

# 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()
Altyapı:

Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)
Altyapı:

new_cases için ham değerleri ve düzleştirilmiş değerleri çizerek düzleştirme sürecinin ne yaptığını görselleştirelim.

🎩 Bazı veri dönüşümleri yaptıktan sonra sonuçları görselleştirmek her zaman iyi bir pratiktir. Bu, sezgi kazanmanızı sağlar ve dönüşümünüzün amaçladığınız şeyi yaptığını akıl süzgecinden geçirmenize yardımcı olur!

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()
Altyapı:

Olabilirlikleri hesaplama

İkinci adım olabilirlikleri hesaplamaktır. Öncekiyle aynı yaklaşımı izleyeceğiz ancak önemli bir farkla: Olabilirlikler yerine log-olabilirlikleri hesaplayacağız. Bu, Kevin'ın önerdiği, $R_t$ hesabında yalnızca son $m$ aralığın kullanılmasına yönelik değişikliği uygulamak için bunları hareketli bir pencerede düzleştirmeyi kolaylaştıracaktır.

Önce olduğu gibi, her gün için olabilirlikleri iç içe bir liste olarak hesaplayacak ve ardından sonuçları tidyr::unnest kullanarak düz bir tabloya açacağız.

🎩 Liste sütunları, verilerle çalışırken son derece kullanışlıdır. Burada kullanılan tasarım deseni pek çok durumda faydalıdır.

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()
Altyapı:

Artgönderimleri hesaplama

Üçüncü adım artgönderim olasılıklarını hesaplamaktır. Hatırlayın

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

Bunu log-olabilirlikler cinsinden şöyle yazabiliriz

$$ 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 paketindeki rollapplyr fonksiyonunu kullanarak log-olabilirliklerin 7 günlük hareketli toplamını hesaplayıp üsselini alarak artgönderimi elde edebiliriz. Son olarak artgönderimleri 1.0'a normalize ederiz.

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()
Altyapı:

Hesapladığımız artgönderim olasılıklarını görselleştirelim. alpha = 0.2 ayarını yaparak üst üste binmeyi azalttığımıza dikkat edin. Bu, kaymakta olan artgönderimleri görselleştirmemize olanak sağlar.

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()
Altyapı:

Rt'yi tahmin et

Son adım, $R_t$ değerlerini ve onları çevreleyen en yüksek yoğunluk aralıklarını tahmin etmektir. HDIInterval::hdi fonksiyonunu uygulamak için artgönderim olasılıklarını kullanarak $r_t$'nin rastgele değerlerini simüle etmemiz gerektiğini unutmayın.

# 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()
Altyapı:

Nihayet gerçeklerle yüzleşme anına geldik! $R_t$'nin tahmini değerlerini görselleştirelim

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()
Altyapı:

Peki bu, Kevin'ın yazısındaki tahminlerle nasıl karşılaştırılıyor? Görünüşe göre modelimdeki ilk tahminler, Kevin'ın tahminlerinden sapıyor. Ancak bu, düzleştirme algoritmalarındaki farkların ya da ilk tahminleri çevreleyen büyük belirsizliğin bir sonucu olabilir. Eğrinin genel şekli ve tepe noktaları oldukça yakın görünüyor. Görev tamam! 🎉

Tüm eyaletlerde döngü

Şimdi tüm eyaletlerde döngü kurup bu tahminleri hesaplama zamanı. Bunu, state alanına göre gruplayarak, veriyi eyalet başına bir tabloya bölüp, her eyalet için $R_t$ tahmini yapıp sonuçları tek bir tabloda birleştirmek üzere purrr::map_df kullanarak kolayca yapabiliriz.

# ⚠️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()
Altyapı:

estimates_all %>%
  display_random()
Altyapı:

Artık tüm eyaletlerdeki tahminlerin küçük çarpanlar grafiğini oluşturabiliriz. ggplot2 kullanımı bunu gerçekten kolaylaştırıyor; tek yapmamız gereken fazladan bir satır kod eklemekti!

# 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)
Altyapı:

Küçük çarpanlar grafiği güzel olsa da, coğrafi konuma göre bir eyalete gitmek zor olabilir. Ya her eyaleti coğrafi konumuna yakın bir yerde çizebilseydik?

Neyse ki, Ryan Hafen tek bir satır kodla bu küçük çarpanları coğrafi bir ızgara üzerine yerleştirmemize olanak tanıyan geofacet paketini oluşturdu.

⚠️ Ne yazık ki, geofacet paketini colab üzerinde kuramadım (sf paketinin gerektirdiği sistem bağımlılıkları nedeniyle olduğunu düşünüyorum); bu nedenle aşağıdaki kod bloğu çalıştırıldığında hata verecektir. Yine de yerel olarak ürettiğim grafiğin kaydedilmiş bir sürümünü ekledim; böylece nasıl göründüğüne dair fikir edinebilirsiniz.

# ⚠️ 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)
Altyapı:

Son olarak, eyaletleri $R_t$'nin en olası tahmini değerine göre sıralayan ve kilitlenme durumuna göre renklendiren Kevin'ın yazısındaki grafiği yeniden oluşturalım.

⚠️ Lütfen bunun tamamen betimleyici bir grafik olduğunu ve kapanmaların etkinliğine dair bir sonuç çıkarılamayacağını unutmayın

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)
Altyapı:

Sonuçlar ve sonraki adımlar

Bu yazıda, Kevin'ın analizini tidyverse kullanarak R'de yeniden ürettik. Olabilirlikleri ve artgönderimleri hesaplamaya yönelik bazı kullanışlı yöntemler öğrendik. Bu süreçte ben de çok şey öğrendim!

Birçok sonraki adım düşünülebilir, ancak aklıma gelen önemli bir tanesini vurgulayacağım. Bu yazıda geliştirilen fonksiyonlar, location, date, cases olmak üzere üç sütunlu oldukça basit bir girdi tablosuna dayandığından, bunu bir R paketine dönüştürmek ilginç bir sonraki adım olabilir. Bu, diğer coğrafyalar için $R_t$ tahminini kolaylaştıracaktır.

Yazımı, ilginç bir yazı yayımladığı ve kodunu paylaştığı için Kevin Systrom'a teşekkür ederek bitirmek istiyorum. Ayrıca bu yazının ilk sürümüne mükemmel geri bildirimler veren DataCamp'teki meslektaşlarıma da teşekkür ederim.

Bu yazı hakkında geri bildiriminiz/yorumlarınız varsa ya da başka genişletmeler önermek isterseniz, bana Twitter'da @ramnath_vaidya ve github üzerinden ulaşabilirsiniz. R programlama hakkında daha fazla bilgi edinmek için