Čekáme na plochou křivku, pojďme si ji zatím nakreslit

Zaslechl jsem zprávu, že cílem karantény a dalších opatření je oploštit křivku. Co tedy dělat, když máme domácí úkoly hotové, látkové roušky ušité a internet pod náporem nadšených homeofficáků kolabuje? Můžeme si kreslit.

Tahle pandemie je skvělá v tom, jak máme takřka real-time data. Nebudu tady hořekovat nad tím, jak jsou nepřesná. Nebo… trochu si tedy zalamentuji:

• každý stát může testovat trochu jinak,
• každý stát zvládá testovat jiné množství lidí,
• některé státy si vyloženě statistiku upravují podle politické situace
• a podobně.

Všude jsou počty nakažených tak trochu podhodnocené, ale v každé zemi v jiné míře. Takže vlastně máme dost přesnou představu pomocí nepřesných čísel. To jen pro začátek, abychom od těchto statistik neměli žádná velká očekávání. Navíc bych chtěl upozornit, že nejsem epidemiolog (což ale lidé v čele této země evidentně také ne).

Nuže k věci. Zahlédl jsem tuhle na Wikipedii krásný obrázek epidemické křivky porovnávající počet nakažených v jednotlivých státech srovnaný ode dne zjištění 100. případu nakažení v dané zemi. Obrázek je bohužel statický a data se každodenně mění, proto by bylo fajn mít něco takového doma a pouštět si to v kormutlivých chvílích těchto dnů.

Získáváme data

Pro tvorbu grafu použijeme R. Aktuální data poskytuje CSSE Johns Hopkins University na GitHubu a jsou dostupná počínaje 22. lednem 2020, kdy už epidemie v Číně běžela naplno. Proto v tomto datasetu Čínu zobrazovat nebudu, v relativním zobrazení by nebyla vypovídající, protože 22. ledna už měla svých 100 případů dávno za sebou (konkrétně 548, jestli se správně dívám). Jsou ovšem i jiné zdroje, jak si ukážeme později.

Začneme tedy načtením balíčků a dat.

library(readr)       # read data
library(dplyr)       # data manipulation
library(tidyr)       # tidy data
library(lubridate)   # parse date
library(ggplot2)     # charts

data <- read_csv("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv")

Upravujeme data

Když se na zdrojová data podíváte, zjistíte, že mají netradiční tvar. Nové dny se nepřipisují dolů (jako položky na účtence), ale doprava. S tím se poněkud hůře pracuje, takže si data převedeme tak, aby byla pod sebou a nikoli vedle sebe. K tomu nám dobře poslouží funkce gather()  z balíčku  tidyr.

Zbavíme se tedy sloupců Province/State, který obsahuje jen části větších států (Austrálie, USA, Čína a podobně), což my nyní rozlišovat nepotřebujeme, a dále Long a Lat, což jsou souřadnice pro případ, že bychom chtěli dělat mapy. To zatím nechceme. V dalším kroku převedu tabulku tak, aby byla data pod sebou (viz výše), přejmenuji sloupec Country/Region na prosté country a sečtu případy pro dané země (to se hodí právě u těch, které mají statistiky rozdělené podle provincií, zkrátka je sečtu do jedné položky za daný stát).

covid_data <- data %>%
  select(-"Province/State", -Long, -Lat) %>%
  gather(date, cases, -"Country/Region") %>%
  rename(country = "Country/Region") %>%
  group_by(country, date) %>%
  summarize(confirmed = sum(cases))

Nyní nám vznikla tabulka se třemi sloupci, tohle jsou první dva řádky:

Všimněte si hrůzostrašného formátu data. Naštěstí je tady funkce mdy() z balíčku  lubridate.

covid_data$date <- covid_data$date %>%
  mdy()

Teď už tabulka vypadá lépe:

Ještě upravím země, které mě budou zajímat. Data z Taiwanu prošla dramatickým vývojem, přechodnou dobu byla zavzata pod Čínu, proto je u Taiwanu hvězdička. A Jižní Korea je zapsána jako „Korea, South“, což můžeme taky změnit. I některé další země by zasloužily upravit pro potřeby grafu, ale ty teď používat nebudu, takže jen tyto dvě. Použijeme k tomu funkci recode() z balíčku  dplyr.

covid_data$country <- covid_data$country %>%
  recode(`Taiwan*` = "Taiwan",
         `Korea, South` = "South Korea")

A už můžeme vyfiltrovat země, které nás zajímají. Je fajn odstartovat zobrazení na 100. případu, protože úvodní fáze nákazy jsou v zemích různě rozpačité. A vyberu si náhodně země, které mne zajímají.

covid_plot <- covid_data %>%
  filter(confirmed > 100 & country %in% c("US", "Taiwan", "Italy", "South Korea", "Czechia", "Singapore", "United Kingdom", "Slovakia", "Spain", "Japan")) %>%
  group_by(country) %>%
  mutate(days100 = as.numeric(date - min(date)))

Chci na konec každé řady dát popisek dané země (místo klasické legendy). Proto si odfiltruji poslední záznamy pro danou zemi a jen ty pak oštítkuji

covid_last <- covid_plot %>%
  slice(which.max(date))

Hurá na graf!

A už můžeme tvořit graf díky balíčku ggplot2. Pro estetický zážitek použiji „vyhlazovací“ funkci geom_smooth(), ale stejně tak by se dala data zobrazit v surovém stavu.

ggplot(covid_plot, aes(days100, confirmed, color = country, label = country)) +
  scale_y_log10(breaks = 10^seq(1, 5)) +
  geom_smooth(se = FALSE) +
  theme_minimal() +
  theme(legend.position = "none") +
  geom_text(data = covid_last,
            hjust = 1, nudge_y = 0.07) +
  labs(x = "Počet dní od 100. případu v dané zemi.", y = NULL,
       title = "Počet nakažených COVID-19 ve vybraných zemích",
       subtitle = paste0("Data ke dni ", format(max(covid_plot$date), "%x"), "."))

Další zdroje dat

Na začátku jsem slíbil, že se podíváme i na jiný zdroj dat. Ten nabízí například balíček nazvaný příznačně nCov19. Je to balíček čínského statistika, takže bohužel tam nenajdete Taiwan. Čísla jsou také lehce odlišná, záleží vždy na zdrojích, odkud autoři čerpají. Balíček v defaultním nastavení stahuje Wuhanský dataset z GitHubu.

Tento kód už nebudu komentovat a vložím jej rovnou s ukázkou grafu:

library(dplyr)       # data manipulation
library(ggplot2)     # charts
library(nCov2019)    # COVID-19 data from https://github.com/canghailan/Wuhan-2019-nCoV

data <- load_nCov2019()

covid_plot <- data["global"] %>%
  as_tibble %>%
  filter(cum_confirm > 100 & country %in% c("China", "Italy", "South Korea", "Czech Republic", "Singapore", "United Kingdom", "Slovakia", "Spain", "Japan")) %>%
  group_by(country) %>%
  mutate(days100 = as.numeric(time - min(time)))

covid_last <- covid_plot %>%
  slice(which.max(time))

ggplot(covid_plot, aes(days100, cum_confirm, color = country, label = country)) +
  scale_y_log10(breaks = 10^seq(1, 4)) +
  geom_smooth(se = FALSE) +
  theme_minimal() +
  theme(legend.position = "none") +
  geom_text(data = covid_last,
            hjust = 0, nudge_x = 0.2) +
  labs(x = "Počet dní od 100. případu v dané zemi.", y = NULL,
       title = "Počet nakažených COVID-19 ve vybraných zemích",
       subtitle = paste0("Data ke dni ", format(time(data), "%x"), "."))

Všimněte si, že tvar křivek je skutečně trochu jiný. Částečně na tom má podíl „vyhlazování“ křivky, ale data u jednotlivých dnů jsou také trochu odlišná, když oba datasety porovnáme. Celkové vyznění grafu ale zůstává stejné.

A je hotovo

Podobně bychom si mohli sestavit relativní grafy (např. počet nakažených v přepočtu na počet obyvatel), posouvat hranici (např. od 10. nebo od 30. případu nakažení) a podobně. Základ kódu ale zůstane víceméně stejný, pokud jsem tedy někde neudělal chybu, což se může stát. :-)

Takže máme nakresleno, můžeme čekat na oploštění křivky!

Aktualizace

• Na podnět RandolfCZ jsem dal zdrojový kód na git a přetvořil skript na Shiny aplikaci. Na radu K> jsem v té aplikaci nepoužil vyhlazování. Díky za tipy!