Wskaźnik zapadalności (dalej WuZet) dla ostatnich 14 dni po naszemu w oryginale zaś 14-day notification rate of newly reported COVID-19 cases per 100 000 population (Data on 14-day notification rate of new COVID-19 cases and deaths). Nasz wspaniały rząd walczący z COVIDem przyjął w swoim czasie `nową strategię' (od tego czasu już kilka razy odnawianą), a w niej był podział na strefy zielona/żółta/czerwona, definiowane odpowiednio jako wartości WuZet poniżej 6 (zielona); 6--12 (żółta) oraz powyżej 12 (czerwona). Dla Sopotu na przykład, który oficjalnie ma około 35 tys mieszkańców, do wejście do czerwonej strefy wystarczały zatem około 4 zakażenia (w ciągu dwóch tygodni). To wszysto wydaje się dziś śmieszne jak ostatnio zakażeń dziennie potrafi być na poziomie trzy-tygodniowej dawki...

Parę dni temu mój bank danych nt. COVID19 uzupełniłem o dane dotyczące liczby zakażeń w powiatach województwa pomorskiego. Akurat WSSE w Gdańsku takie dane w sposób w miarę porządny publikuje i da się je względnie łatwo odzyskać z raportów publikowanych i archiwizowanych (brawo--na to nie wpadł nawet Minister w MZ) pod adresem http://www.wsse.gda.pl/.

library("dplyr")
library("tidyr")
library("ggplot2")
m1unit <- 100000
# pomorskie_by_powiat.csv pobrane z www.wsse.gda.pl
# format: data;powiat;nc 
d <- read.csv("pomorskie_by_powiat.csv", sep = ';',  header=T, na.string="NA" )

# wartości skumulowane (po powiatach dlatego tak dziwnie)
# replace_na jest potrzebne bo cumsum nie obsługuje NA
e <- d %>% group_by(powiat) %>% dplyr::mutate(tc = cumsum(replace_na(nc, 0)),
tc1m = cumsum(replace_na(nc, 0)) / pop * m1unit
)
## dzień ostatni
day00 <- as.Date(last.obs)
## dwa tygodnie przed ostatnim
day14 <- as.Date(last.obs) - 14
## 4 tygodnie przed ostatnim
day28 <- as.Date(last.obs) - 28

## Stan na dzień ostatni
e00 <- e %>% filter (as.Date(dat) == day00 ) %>%
 group_by(powiat) %>% as.data.frame

## BTW Żeby było dziwniej zapis
## e0 <- e %>% group_by(powiat) %>%  slice_tail(n=1) %>% as.data.frame
## Daje inny wynik, liczby te same, ale porządek wierszy jest inny

## Stan na dzień dwa tygodnie przed
e14 <- e %>% filter (as.Date(dat) == day14 ) %>%
 group_by(powiat) %>% as.data.frame

e14 <- e %>% filter (as.Date(dat) == day14 ) %>%
 group_by(powiat) %>% as.data.frame

## Stan na dzień 4 tygodnie przed
e28 <- e %>% filter (as.Date(dat) == day28 ) %>%
 group_by(powiat) %>% as.data.frame

## nowe zakażenia w tygodniach 3/4
c28 <- e14$tc - e28$tc

## nowe zakażenie w tygodniach 2/1
c14 <- e00$tc - e14$tc
## To samo co c14/c28 ale w przeliczeniu
## na 100 tys:
c28m1 <- e14$tc1m - e28$tc1m
c14m1 <- e00$tc1m - e28$tc1m
## Dynamika zmiana WuZet w dwóch ostatnich okresach
## tj WuZet12 względem WuZet34
#d14v28 <- (c14m1 - c28m1) / c28m1 * 100
d14v28 <- (c14m1/c28m1) * 100
##
## Można sobie teraz c14m1/d14v28 na wykresach przestawić

Na dzień 7 listopada wyniki były takie (jeżeli się nie rąbnąłem w powyższym kodzie):

Na dzień 7 listopada wyniki były takie (jeżeli się nie rąbnąłem w powyższym kodzie):

sprintf("%10.10s = %6.2f | %6.2f | %6.2f - %6.2f | %4i | %4i | %4i (%4i)",
   e00$powiat, d14v28, c14m1, e00$tc1m, e28$tc1m, e00$tc, e14$tc, e28$tc, e00$pop )
 [1] "    Gdynia = 229.15 | 577.74 | 1129.08 - 299.22 | 2781 | 1358 |  737 (246306)"
 [2] "    Gdańsk = 149.50 | 416.37 | 1045.98 - 351.10 | 4856 | 2923 | 1630 (464254)"
 [3] "    Słupsk = 228.26 | 803.59 | 1387.42 - 231.78 | 1269 |  534 |  212 (91465)"
 [4] "     Sopot = 144.90 | 583.03 | 1404.21 - 418.80 |  513 |  300 |  153 (36533)"
 [5] "  tczewski = 437.50 | 905.98 | 1323.60 - 210.53 | 1534 |  484 |  244 (115896)"
 [6] "   gdański = 399.21 | 889.61 | 1353.71 - 241.26 | 1543 |  529 |  275 (113983)"
 [7] "  kartuski = 197.07 | 855.49 | 1846.22 - 556.63 | 2471 | 1326 |  745 (133841)"
 [8] "  bytowski = 268.71 | 998.31 | 1693.33 - 323.50 | 1340 |  550 |  256 (79134)"
 [9] " malborski = 329.61 | 923.16 | 1408.21 - 204.97 |  900 |  310 |  131 (63911)"
[10] "     pucki = 225.35 | 766.35 | 1545.69 - 439.26 | 1309 |  660 |  372 (84687)"
[11] "wejherowsk = 150.71 | 396.16 |  971.92 - 312.90 | 2078 | 1231 |  669 (213803)"
[12] "starogardz = 216.36 | 744.52 | 1388.93 - 300.31 | 1776 |  824 |  384 (127868)"
[13] " chojnicki = 266.33 | 813.36 | 1311.04 - 192.29 | 1275 |  484 |  187 (97251)"
[14] "  sztumski = 309.52 | 619.84 |  979.83 - 159.73 |  411 |  151 |   67 (41946)"
[15] " kwidzyńsk = 251.34 | 563.39 |  973.35 - 185.80 |  812 |  342 |  155 (83423)"
[16] " kościersk = 293.30 | 786.88 | 1392.60 - 337.43 | 1007 |  438 |  244 (72311)"
[17] "nowodworsk = 263.21 | 777.35 | 1273.30 - 200.61 |  457 |  178 |   72 (35891)"
[18] "   słupski = 244.23 | 514.50 |  969.24 - 244.08 |  957 |  449 |  241 (98737)"
[19] "  lęborski = 172.27 | 618.09 | 1135.18 - 158.29 |  753 |  343 |  105 (66333)"
[20] " człuchows = 268.29 | 388.16 |  571.65 -  38.82 |  324 |  104 |   22 (56678)"

Dynamika WuZeta: Gdynia/Gdańsk/Sopot = 129.1%, 149.5% i 144.9%; wartość maksymalna 437% (tczewski); wartość minimalna 150% (wejherowski). Najnowsze wartości współczynnika WuZet: Gdynia/Gdańsk/Sopot= 577.7, 416.4 oraz 583.0; wartość maksymalna 998,3 (bytowski); wartość minimalna 388.1 (człuchowski). Dane i skrypty są tutaj.

Skrypt rysujący wykres łącznej liczby zakażeń i zgonów wg kontynentów. Przy okazji podjąłem nieśmiajłą próbę zapanowania nad (jednolitym) wyglądem. Udając, że istnieje coś takiego jak Narodowy Instytut Korona Wirusa stworzyłem dla niego wizualizację, którą dodaje się do wykresu w postaci jednego polecenia. To jedno polecenie to funkcja theme_nikw:

#!/usr/bin/env Rscript
#
library("ggplot2")
library("dplyr")
library("scales")
library("ggthemes")
library("ggpubr")
#
theme_nikw <- function(){
 theme_wsj() %+replace%
 theme(
  panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "cornsilk3"),
  panel.grid.major = element_line(size = 0.25, linetype = 'solid', colour = "cornsilk3"),
  ##
  axis.text.x  = element_text(size = 6 ),
  axis.text.y  = element_text(size = 6 ),
  ## https://stackoverflow.com/questions/14379737/how-can-i-make-xlab-and-ylab-visible-when-using-theme-wsj-ggthemes
  axis.title  = element_text(family="sans", size=6)
  ## Poniższe natomiast nie działa:
  #axis.title.x  = element_text(size = 6 ),
  #axis.title.y  = element_text(size = 6 ),
  ## margin=margin(r=, t=, b = 3, l=, unit="pt")
  plot.title=element_text(family="sans", size=14, hjust=0, margin=margin(b = 3, unit="pt")),
  plot.subtitle=element_text(family="sans", size=8, hjust=0),
  legend.title=element_text(family="sans", size=8),
  plot.caption = element_text(family="sans", size = 6)
  )
}

Uprzedzając wydarzenia, kolejna funkcja notin służy do tego co funkcja in tyle, że zamiast zwracać wartości pasujące, zwraca te które nie pasują:

## https://stackoverflow.com/questions/34444295/how-to-specify-does-not-contain-in-dplyr-filter
`%notin%` = function(x,y) !(x %in% y)

Teraz ustawiam różne wartości globalne (options służy do wyłączenia zapisu liczb w trybie scientic notation).

options(scipen = 999) ## albo options(scipen = 999, digits=3)
note <- "(c) NI-KW @ github.com/knsm-psw/NI-KW"
today <- Sys.Date()
tt <- format(today, "%Y-%m-%d")
spanV <- 0.25
mainBreaks <- "4 weeks"
mainColor <- "deeppink"
loessColor <- "deeppink"

## Przed/po \n musi być odstęp inaczej nie działa
## url <- "https://www.ecdc.europa.eu/en/publications-data/ \n download-todays-data-geographic-distribution-covid-19-cases-worldwide"
url <- "www.ecdc.europa.eu/en/publications-data/download-todays-data-geographic-distribution-covid-19-cases-worldwide"
surl <- sprintf ("Data retrived from %s", url)

Czytanie zbiorów danych. Dane dotyczące zakażeń pochodzą ze strony ECDC (European Center for Disease Prevention and Control). Dane dotyczące ludności też z tej strony tyle, że wydłubałem je raz i trzymam w oddzielnym pliku zamiast powielać niepotrzebnie w pliku z danymi o COVID19 (jak to robi ECDC)

## info o wielkości populacji i kontynencie
c <- read.csv("ecdc_countries.csv", sep = ';',  header=T, na.string="NA" )
## dane dzienne/cały świat (ECDC)
d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", 
   colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));

Złączenie zbiorów w oparciu o kolumnę id:

d <- left_join(d, c, by = "id")

Parę rzeczy trzeba uporządkować. Po pierwsze usunąć wpisy zawierające kraj bez ISO-kodu (coś jest nie tak z transformacją XSLX na CSV, muszę to sprawdzić). Po drugie czasami liczby zakażeń/śmierci są ujemne bo dane zostają zrewidowane. Nie jest to najlepsze rozwiązanie, ale żeby przynajmniej na pierwszy rzut oka wykres nie wyglądał absurdalnie to zamieniam liczby ujemne na NA:

# czyszczenie
# tylko kraje z ID
d <- d %>% filter(!is.na(id))
# zamień na liczby (był problem z czytaniem CSV bezpośrednio)
d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)

# Dane bezsensowne zamień na NA
# Zdarzają się ujemne bo kraje zmieniają klasyfikację/przeliczają
d$newc[ (d$newc < 0) ] <- NA
d$newd[ (d$newd < 0) ] <- NA

Liczę wartości na 1/mln. Obliczam dzienne wartości zakażeń/śmierci na 1mln wg kontynentów

## zgony na 1mln
d$tot1m <- d$totald/d$popData2019 * 1000000
## Oblicz sumy dla kontynentów
## Kontynentów jest więcej niż 5:-) w zbiorze (jakieś śmieci)
continents <- c('Africa', 'America', 'Asia', 'Europe', 'Oceania')
## ogółem wg 
cont <- d %>% filter (continent %in% continents) %>% group_by(date,continent) %>% summarise(
  tc = sum(newc, na.rm=TRUE),
  td = sum(newd, na.rm=TRUE),
  tc1m = sum(newc, na.rm=TRUE)/sum(popData2019, na.rm=TRUE) * 1000000,
  td1m = sum(newd, na.rm=TRUE)/sum(popData2019, na.rm=TRUE) * 1000000
  ) %>% as.data.frame
## str(tc)

## Z ciekawości
noncont <- d %>% filter (continent %notin% continents) %>% as.data.frame
print(noncont)

Wykreślam wykresy. Punkty przestawiają dane dzienne. Krzywa loess trend. Źródło danych jest w podtytule. W napisie umieszczanym pod wykresem jest informacja o autorze i link do repozytorium githuba (polecenie labs(caption=note))

pd1 <- ggplot(cont, aes(x= as.Date(date, format="%Y-%m-%d"), color = continent, y=tc)) +
 geom_point(aes(y=tc, color = continent), size=.4, alpha=.5) +
 geom_smooth(method="loess", se=F, span=spanV, aes(fill=continent, y=tc, group=continent)) +
 xlab(label="") + ylab(label="cases") + 
 scale_x_date( labels = date_format("%m/%d"), breaks = mainBreaks) +
 theme_nikw() +
 labs(caption=note) +
 ggtitle(sprintf ("COVID19: cases (%s)", last.obs), subtitle=sprintf("%s", surl))

pd2 <- ggplot(cont, aes(x= as.Date(date, format="%Y-%m-%d"), , color = continent, y=td)) +
 geom_point(aes(y=td, color = continent), size=.4, alpha=.5) +
 geom_smooth(method="loess", se=F, span=spanV, aes(fill=continent, y=td, group=continent)) +
 xlab(label="") + ylab(label="deaths") + 
 scale_x_date( labels = date_format("%m/%d"), breaks = mainBreaks) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 theme_nikw() +
 labs(caption=note) +
 ggtitle(sprintf ("COVID19: deaths (%s)", last.obs), subtitle=sprintf("%s", surl))

ggsave(plot=pd1, file="tc_by_c.png", width=10)
ggsave(plot=pd2, file="td_by_c.png", width=10)

pd1m <- ggplot(cont, aes(x= as.Date(date, format="%Y-%m-%d"), color = continent, y=tc1m)) +
 geom_point(aes(y=tc1m, color = continent), size=.4, alpha=.5) +
 geom_smooth(method="loess", se=F, span=spanV, aes(y=tc1m )) +
 xlab(label="") + ylab(label="cases") + 
 scale_x_date( labels = date_format("%m/%d"), breaks = mainBreaks) +
 theme_nikw() +
 labs(caption=note) +
 ggtitle(sprintf ("COVID19: cases per 1mln (%s)", last.obs), subtitle=sprintf("%s", surl))

pd2m <- ggplot(cont, aes(x= as.Date(date, format="%Y-%m-%d"), color = continent, y=td1m)) +
 geom_point(aes(y=td1m, color = continent), size=.4, alpha=.5) +
 geom_smooth(method="loess", se=F, span=spanV, aes(y=td1m )) +
 xlab(label="") + ylab(label="deaths") + 
 scale_x_date( labels = date_format("%m/%d"), breaks = mainBreaks) +
 theme_nikw() +
 labs(caption=note) +
 ggtitle(sprintf ("COVID19: deaths per 1mln(%s)", last.obs), subtitle=sprintf("%s", surl))

ggsave(plot=pd1m, file="tc1m_by_c.png", width=10)
ggsave(plot=pd2m, file="td1m_by_c.png", width=10)

Wszystko. Skrypt jest wykonywany automatycznie o ustalonej godzinie, po czym wykresy wysyłane są na twitterowe konto "Instytutu".

Dokument #1

Dokument #2

Pobielenie

Dokument gotowy do OCR

Że się pandemia rozwija zachciało mi się odzyskać dane podawane przez polskie Ministerstwo Zdrowia na Twitterze w formie komunikatów obrazkowych. Otóż raz dziennie o 10:00 pojawia się rysunek zawierający dane dotyczące zajętych łóżek dedykowanych COVID-19, zajętych respiratorów, liczby osób objętych kwarantanną i jeszcze dwóch innych rzeczy. Ponadto, raz na tydzień w poniedziałek pojawia się rysunek z liczbą testów w podziale na województwa. Tych danych nie ma na stronie MZ, a jeżeli są to tak publikowane, że ja nie potrafię tego odszukać. Że MZ udostępnia dane obrazkowo za pośrednictwem amerykańskiej firmy zamiast w sposób opisany ustawą o dostępnie do informacji publicznej to oczywiście skandal i żenada.

Cały stream https://twitter.com/MZ_GOV_PL to ja mam ściągnięty. Tam nie ma obrazków; są URLe do obrazków, które można łatwo pobrać. Ja to robię na raty, najpierw json do csv:

#!/usr/bin/perl
# Zamiana Json -> CSV
use JSON;
use Time::Piece;
use open ":encoding(utf8)";
use open IN => ":encoding(utf8)", OUT => ":utf8";
binmode(STDOUT, ":utf8");

print "id;date;repid;text\n";

while (<>) { chomp();
  $tweet =  $_;

  my $json = decode_json( $tweet );
  $tid = $json->{"id"};
  $dat = $json->{"created_at"};
  $dat = Time::Piece->strptime($dat,
   "%a %b %d %H:%M:%S %z %Y")->strftime("%Y-%m-%dT%H:%M:%S");

  $mmm = $json->{"entities"}{"media"}; ## lista-haszy

  for $mm ( @{$mmm} ) { $media = $mm->{media_url} ;
    ## id-tweeta;data;url-do-rysunku
    print "$tid;$dat;$media\n";  }
}

Teraz można ściągnąć rysunku zwykłym wgetem albo curlem, pobierając URLa z cvsa. Rysunki o łóżkach i wentylatorach pojawiają się codziennie około 10:30 (czyli 8:30 GMT). Rysunki o liczbie testów wg województw w poniedziałki generalnie około 16:00. Więc prostym skryptem Perla ściągam poniedziałkowe rysunki opublikowane po 13:30GMT oraz wszystkie po 8:00GMT a przed 10:00GMT. Po ściągnięciu oglądam i wywalam nierelewantne.

OCR robię programem tesseract:

tesseract PLIK.png PLIK

Powstanie PLIK.txt a w nim tekst z rysunku. Z danymi poniedziałkowymi był problem, tesseract się gubił i PLIK.txt nic nie zawierał Żeby mu pomóc najpierw upraszczałem rysunek:

#!/bin/bash
## pomaluj na biało fragmenty nie zawierające liczb
convert "$1" -fill white -draw "rectangle 0,0 1200,115" \
  -draw "rectangle 0,640 200,675"  \
  -draw "rectangle 350,115 830,540" PLIK_0.png
## zamień wszystkie kolory na biały za wyjątkiem czarnego:
convert PLIK_0.png -fuzz 30% -fill white +opaque black PLIK_1.png
## zrób OCR
tesseract PLIK_1.png PLIK_1
## oczyść i dopisz wiersz do WYNIKI.txt
grep '[0-9]' PLIK_1.txt | \
awk '{gsub(/[ \t]/, ""); l = l ";" $0 }; END{print l}' > WYNIKI.txt

Jeżeli powyższy skrypt nazywa się png2txt.sh to teraz:

for x in *.png; do
   if [ -f $x ] ; then
        png2txt.sh $x
   fi
done

Raporty poniedziałkowe zaczęły być wysyłane od 11 maja 2020. Raporty codzienne pobrałem od początku lipca. W szczególności rysunki `poniedziałkowe' zawierają dane kumulowane. Kiedy na podstawie tych danych utworzyłem dane `tygodniowe' (jako różnica między stanem na bieżący tydzień minus stanem na poprzedni tydzień), to dla województwa świętokrzyskiego i raportu z 10 sierpnia rezultat okazał się ujemny i do tego ogromny. Licznik cofnęło mówiąc kolokwialnie.

Po wpisaniu do google testy+świętokrzyskie się okazało, że sprawa jest znana: przez 2 miesiące województwo świętokrzyskie podawało dane ewidentnie z sufitu raportując 100% wzrost tydzień/tydzień przy maksymalnym dla następnego województwa poziomie 15%... Jak po 2 miesiącach tej twórczości doszli do absurdalnej liczby ktoś się w MZ połapał i napisał (na Twitterze), że należy odjąć te lipne 230 tys...

No to tak z grubsza wygląda rzetelność danych n/t COVID19 w PL...

Pozyskane dane są: tutaj.

W mediach podają albo to co było wczoraj, albo razem od początku świata. Zwłaszcza to łącznie od początku ma niewielką wartość, bo są kraje gdzie pandemia wygasa a są gdzie wręcz przeciwnie. Co z tego, że we Włoszech liczba ofiar to 35 tysięcy jak od miesięcy jest to 5--25 osób dziennie czyli tyle co PL, gdzie zmarło do tej pory kilkanaście razy mniej chorych. Podobnie w Chinach, gdzie cała afera się zaczęła, zmarło we wrześniu 15 osób. Media ponadto lubią duże liczby więc zwykle podają zmarłych, a nie zmarłych na 1mln ludności. Nawet w danych ECDC nie ma tego wskaźnika zresztą. Nie ma też wskaźnika zgony/przypadki, który mocno niefachowo i zgrubie pokazywałby skalę zagrożenia (im większa wartość tym więcej zarażonych umiera)

Taka filozofia stoi za skryptem załączonym niżej. Liczy on zmarłych/1mln ludności oraz współczynnik zgony/przypadki dla 14 ostatnich dni. Wartości przedstawia na wykresach: punktowym, histogramie oraz pudełkowym. W podziale na kontynenty i/lub kraje OECD/pozostałe. Ponieważ krajów na świecie jest ponad 200, to nie pokazuje danych dla krajów w których liczba zmarłych jest mniejsza od 100, co eliminuje małe kraje albo takie, w których liczba ofiar COVID19 też jest mała.

Wynik obok. Skrypt poniżej

#!/usr/bin/env Rscript
#
library("ggplot2")
library("dplyr")
#

## Space before/after \n otherwise \n is ignored
url <- "https://www.ecdc.europa.eu/en/publications-data/ \n download-todays-data-geographic-distribution-covid-19-cases-worldwide"
surl <- sprintf ("retrived %s from %s", tt, url)

## population/continent data:
c <- read.csv("ecdc_countries_names.csv", sep = ';',  header=T, na.string="NA" )

## COVID data
d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", 
   colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));

# today <- Sys.Date()
# rather newest date in data set:
today <- as.Date(last(d$date))
# 14 days ago
first.day <- today - 14
tt<- format(today, "%Y-%m-%d")
fd<- format(first.day, "%Y-%m-%d")

# select last row for every country
##t <- d %>% group_by(id) %>%  top_n(1, date) %>% as.data.frame
## top_n jest obsolete w nowym dplyr
t <- d %>% group_by(id) %>%  slice_tail(n=1) %>% as.data.frame

# cont. select id/totald columns, rename totald to alld
t <- t %>% select(id, totald) %>% rename(alld = totald) %>% as.data.frame
# join t to d:
d <- left_join(d, t, by = "id")
## extra column alld = number of deaths at `today'
##str(d)

# only countries with > 99 deaths:
d <- d %>% filter (alld > 99 ) %>% as.data.frame

## OECD membership data
o <- read.csv("OECD-list.csv", sep = ';',  header=T, na.string="NA" )

## Remove rows with NA in ID column (just to be sure)
## d <- d %>% filter(!is.na(id))

d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)

cat ("### Cleaning newc/newd: assign NA to negatives...\n")
d$newc[ (d$newc < 0) ] <- NA
d$newd[ (d$newd < 0) ] <- NA

## Filter last 14 days only 
d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > fd) %>% as.data.frame
## Last day
last.obs <- last(d$date)
## Format lable firstDay--lastDay (for title):
last14 <- sprintf ("%s--%s", fd, last.obs)

## Sum-up cases and deaths
t <- d %>% group_by(id) %>% summarise( tc = sum(newc, na.rm=TRUE), td = sum(newd, na.rm=TRUE)) %>% as.data.frame 
## Add country populations
t <- left_join(t, c, by = "id")
## Add OECD membership
## if membership column is not NA = OECD member
t <- left_join(t, o, by = "id") 

t$access[ (t$access > 0) ] <- 'oecd'
t$access[ (is.na(t$access)) ] <- 'non-oecd'
str(t)

## Deaths/Cases ratio %%
t$totr <- t$td/t$tc * 100
## Deaths/1mln
t$tot1m <- t$td/t$popData2019 * 1000000

## PL row
dPL <- t[t$id=="PL",]
##str(dPL)
## Extract ratios for PL
totr.PL <- dPL$totr
totr.PL
tot1m.PL <- dPL$tot1m

## Set scales
pScale <- seq(0,max(t$totr, na.rm=T), by=1)
mScale <- seq(0,max(t$tot1m, na.rm=T), by=10)

## Graphs
## 1. deaths/cases ratio (dot-plot)
## color=as.factor(access)): draw OECD/nonOECD with separate colors
p1 <- ggplot(t, aes(x = reorder(country, totr) )) +
  ### One group/one color:
  ##geom_point(aes(y = totr, color='navyblue'), size=1) +
  ### Groups with different colors:
  geom_point(aes(y = totr, color=as.factor(access)), size=1) +
  xlab("country") +
  ylab("%") +
  ggtitle(sprintf("Covid19 death cases ratio (%s)", last14),
    subtitle="Only countries with 100 deaths and more | pale red line = PL level") +
  theme(axis.text = element_text(size = 4)) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ##theme(legend.position="none") +
  scale_color_discrete(name = "Membership") +
  ## Indicate PL-level of  d/c ratio:
  geom_hline(yintercept = totr.PL, color="red", alpha=.25, size=1) +
  labs(caption=surl) +
  scale_y_continuous(breaks=pScale) +
  ##coord_flip(ylim = c(0, 15))
  coord_flip()

ggsave(plot=p1, file="totr_p14.png", width=10)

## 2. deaths/cases ratio (histogram)
p2 <- ggplot(t, aes(x = totr) ) + 
 geom_histogram(binwidth = 0.5, fill='navyblue', alpha=.5) +
 ###
 ylab("N") +
 xlab("%") +
 ggtitle(sprintf("Covid19 deaths/cases ratio (%s)", last14),
   subtitle="Only countries with 100 deaths and more | pale red line = PL level") +
 scale_x_continuous(breaks=pScale) +
 scale_y_continuous(breaks=seq(0, 40, by=2)) +
 geom_vline(xintercept = totr.PL, color="red", alpha=.25, size=1) +
 ##coord_cartesian(ylim = c(0, 30), xlim=c(0, 8))
 labs(caption=surl)

ggsave(plot=p2, file="totr_h14.png", width=10)

## 3. deaths/1m (dot-plot)
## color=as.factor(continent)): draw continents with separate colors
p3 <- ggplot(t, aes(x = reorder(country, tot1m) )) +
  geom_point(aes(y = tot1m, color=as.factor(continent)), size =1 ) +
  xlab("") +
  ylab("deaths") +
  ggtitle(sprintf("Covid19 deaths per 1 million (%s)", last14),
    subtitle="Only countries with 100 deaths and more | red pale line = PL level") +
  theme(axis.text = element_text(size = 6)) +
  geom_hline(yintercept = tot1m.PL, color="red", alpha=.25, size=1) +
  labs(caption=surl) +
  scale_color_discrete(name = "Continent") +
  scale_y_continuous(breaks=mScale) +
  coord_flip()

ggsave(plot=p3, file="totr_m14.png", height=10)

## 3. deaths/1m (box-plots for continents)
p4 <- ggplot(t, aes(x=as.factor(continent), y=tot1m, fill=as.factor(continent))) + 
 geom_boxplot() +
 ylab("deaths") +
 xlab("continent") +
 ggtitle(sprintf("Covid19 deaths per 1 million (%s)", last14),
  subtitle="Only countries with 100 deaths and more") +
 labs(caption=surl) +
 scale_y_continuous(breaks=mScale) +
 theme(axis.text = element_text(size = 6)) +
 theme(legend.position="none")

ggsave(plot=p4, file="totr_c14.png", width=10)

## Mean/Median/IQR per continent
ts <-  t %>% group_by(as.factor(continent)) %>% 
    summarise(Mean=mean(tot1m), Median=median(tot1m), Iqr=IQR(tot1m))
ts

Wartość zmarli/przypadki zależy ewidentnie od słynnej liczby testów. Taki Jemen wygląda, że wcale nie testuje więc ma mało przypadków. Z drugiej strony ci co dużo testują mają dużo przypadków, w tym słynne bezobjawowe więc też nie bardzo wiadomo co to oznacza w praktyce. Bo chyba się testuje po coś, a nie żeby testować dla samego testowania. Więc: dużo testów -- na mój niefachowy rozum -- to wczesne wykrycie, a mało to (w teorii przynajmniej) dużo nieujawnionych/nieizolowanych zarażonych, co powinno dawać w konsekwencji coraz większą/lawinowo większą liczbę przypadków, a co za tym idzie zgonów a tak nie jest.

Konkretyzując rysunek: Ameryka na czele liczby zgonów z medianą prawie 30/1mln, Europa poniżej 3,5/1mln a Azja coś koło 1,4/1mln (średnie odpowiednio 26,4/6,00/5,1.)

Przy okazji się okazało, że mam dplyra w wersji 0.8 co jest najnowszą wersją w moim Debianie. No i w tej wersji 0.8 nie ma funkcji slice_head() oraz slice_tail(), która jest mi teraz potrzebna. Z instalowaniem pakietów R w Linuksie to różnie bywa, ale znalazłem na stronie Faster package installs on Linux with Package Manager beta release coś takiego:

RStudio Package Manager 1.0.12 introduces beta support for the binary format of R packages. In this binary format, the repository's packages are already compiled and can be `installed' much faster -- `installing' the package amounts to just downloading the binary!

In a configured environment, users can access these packages using `install.packages' without any changes to their workflow. For example, if the user is using R 3.5 on Ubuntu 16 (Xenial):

> install.packages("dplyr", repos = "https://demo.rstudiopm.com/cran/__linux__/xenial/latest")

Zadziałało, ale (bo jest zawsze jakieś ale) powtórzenie powyższego na raspberry pi nie zadziała. Wszystko się instaluje tyle że architektury się nie zgadzają. Więc i tak musiałem wykonać

> install.packages("dplyr")
## sprawdzenie wersji pakietów
> sessionInfo()

Instalacja trwała 4 minuty i zakończyła się sukcesem.

Jeszcze jeden skrypt do rysowania danych nt COVID19:

#!/usr/bin/env Rscript
# Przekazywanie argumentów z wiersza poleceń
# np.: Rscript --vanilla c19.R --iso PL -clean
library("optparse")
#
library("ggplot2")
library("dplyr")
library("scales")
library("ggpubr")
#
# parametr wygładzania (loess)
spanV <- 0.25
# UWAGA: przed/po \n musi być odstęp inaczej nie działa
surl <- "https://www.ecdc.europa.eu/en/publications-data/ \n download-todays-data-geographic-distribution-covid-19-cases-worldwide"
c0 <- 'PL'

option_list <- list(
  make_option(c("-i", "--iso"), action="store", type="character", default=c0, help="country ISO2 code"),
  make_option(c("-c", "--clean"), action="store_true", default=T, help="extra clean data?")
); 
 
opt_parser <- OptionParser(option_list=option_list);
opt <- parse_args(opt_parser);

c0 <- opt$iso
dataClean <- opt$clean

# wczytanie danych
# date;id;country;newc;newd;totalc;totald
d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", 
   colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));

str(d)

d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)

## Dane zawierają wartości ujemne co jest bez sensu
## z opcją --clean te ujemne wartości są zamieniane na NA
if ( dataClean ) {
  cat ("### Cleaning newc/newd: assign NA to negatives...\n")
  d$newc[ (d$newc < 0) ] <- NA
  d$newd[ (d$newd < 0) ] <- NA
}

## Współczynniki zmarli/zakażeni
d$newr <- d$newd/d$newc * 100
d$totr <- d$totald/d$totalc * 100

## Wartości > 50% są zamieniane na NA (zwykle >50% wynika z błędnych danych)
if ( dataClean ) {
  cat ("### Cleaning newc/newd: assign NA to newr/totr higher than 50...\n")
  d$newr[ (d$newr > 50) ] <- NA
  d$totr[ (d$totr > 50) ] <- NA
}

## Pomiń obserwacje wcześniejsze niż 15/02
d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame

d0 <- d %>% filter (id == c0) %>% as.data.frame
t0 <- d0 %>% group_by(id) %>%  summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T))

lab0c <- toString(paste (sep=" = ", t0$id, t0$cc))
lab0d <- toString(paste (sep=" = ", t0$id, t0$dd))

## koniecznie dodać na.rm=T bo inaczej zwraca NA (jeżeli znajdzie NA)
maxCC <- max (d0$newc, na.rm=T)
maxDD <- max (d0$newd, na.rm=T)
maxRR <- max (d0$totr, na.rm=T)

last.obs <- last(d0$date)
first.date <- first(d0$date)
fstDay <- as.Date(first.date)
last.totr <- last(d0$totr)
max.newr <- max(d0$newr, na.rm=T)

## Przykład dodania 14 dni do daty
## srcDay <- as.Date(first.date) +14
## https://stackoverflow.com/questions/10322035/r-adding-days-to-a-date
srcDay <- as.Date(last.obs)

## Nazwa pliku wynikowego
## c19_ISO2_DATA.png, gdzie DATA jest datą ostatniej obserwacji
## np.: c19_SE_2020-09-16.png
c0o <- sprintf ("c19_%s_%s.png", c0, last.obs)

## Rysunek1: nowe przypadki
pc0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=newc)) + 
 geom_point(aes(group = id, color = id), size=.8) +
 geom_smooth(method="loess", se=F, span=spanV) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 annotate("text", x = fstDay, y = 0.95 * maxCC, 
  label = sprintf("Total: %i cases", t0$cc), hjust = 0, vjust=0,
  alpha=0.3, color='steelblue', size=6) +
  xlab(label="") +
 ## Nie drukuj legendy
 theme(legend.position="none") +
 ggtitle(sprintf("%s: new confirmed cases (%s)", t0$id, last.obs))

## Rysunek2: nowe zgony
pd0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=newd)) + 
 geom_point(aes(group = id, color = id), size=.8) +
 geom_smooth(method="loess", se=F, span=spanV) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 annotate("text", x = fstDay, y = 0.95 * maxDD, label = sprintf("Total: %i deaths", t0$dd), hjust = 0, vjust=0,
  alpha=0.3, color='steelblue', size=6) +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 xlab(label="") +
 theme(legend.position="none") +
 ggtitle(sprintf ("%s: deaths (%s)", t0$id, last.obs))

## Rysunek3: nowe zgony/przypadki *100%
pr0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=newr)) + 
 geom_point(aes(group = id, color = id), size=.8) +
 geom_smooth(method="loess", se=F, span=spanV) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 annotate("text", x = fstDay, y = 0.95 * max.newr, label = sprintf("Maximum: %.2f %%", max.newr), hjust = 0, vjust=0,
  alpha=0.3, color='steelblue', size=6) +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 xlab(label="") +
 ylab(label="%") +
 theme(legend.position="none") +
 ggtitle(sprintf ("%s: deaths/cases %% (%s)", t0$id, last.obs) )

## Rysunek4: łączne zgony/przypadki *100%
prt0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=totr)) + 
 geom_point(aes(group = id, color = id), size=.8) +
 geom_smooth(method="loess", se=F, span=spanV) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 annotate("text", x = fstDay, y = 0.95 * maxRR, label = sprintf("Average: %.2f %%", last.totr), hjust = 0, vjust=0,
  alpha=0.3, color='steelblue', size=6) +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 xlab(label="") +
 ylab(label="%") +
 theme(legend.position="none") +
 annotate("text", x = srcDay, y = 0, label = surl, hjust = 1, alpha=.3, size=3) +
 ggtitle(sprintf ("%s: total deaths/cases %% (%s)", t0$id, last.obs) )

p00 <- ggarrange(pc0,pd0, pr0, prt0, ncol=2, nrow=2)

ggsave(plot=p00, c0o, width=15)

Użycie:

Rscript --vanilla c19.R --iso PL

albo:

for i in 'AU' 'BR' 'IN' 'US' 'ES' 'SE' 'PL' 'DE' 'GB'; 
 Rscript --vanilla c19.R --iso $i
done

Przykładowy wynik dla Hiszpanii

Jak widać jakość danych jest katastrofalna: pojawiające się liczne zera to w rzeczywistości brak danych. Zwykle sobota/niedziela zero a potem sruuuu 30 tysięcy bo za trzy dni razem. Wszyscy są zmęczeni w tej Hiszpanii pandemię i nawet nie chce im się danych podsyłać do ECDC?

MZ to Ministerstwo Zdrowia. Specyfiką PL jest brak danych publicznych nt. Pandemii.. Na stronie GIS nie ma nic, a stacje wojewódzkie publikują jak chcą i co chcą. Ministerstwo zdrowia z kolei na swojej stronie podaje tylko dane na bieżący dzień, zaś na amerykańskim Twitterze publikuje komunikaty na temat. To co jest na stronie znika oczywiście następnego dnia zastępowane przez inne cyferki. Do tego są tam tylko dzienne dane zbiorcze o liczbie zarażonych i zmarłych, (w podziale na województwa). Nie ma na przykład wieku, płci itp... To co jest na Twitterze z kolei ma formę tekstowego komunikatu postaci: Mamy 502 nowe i potwierdzone przypadki zakażenia #koronawirus z województw: małopolskiego (79), śląskiego (77), mazowieckiego (75)... [...] Z przykrością informujemy o śmierci 6 osób zakażonych koronawirusem (wiek-płeć, miejsce zgonu): 73-K Kędzierzyn-Koźle, 75-M Łańcut, 92-K Lipie, 72-M, 87-M, 85-K Gdańsk.

Czyli podają zarażonych ogółem i w podziale na województwa oraz dane indywidualne zmarłych w postaci płci i wieku oraz miejsca zgonu (miasta żeby było inaczej niż w przypadku podawania zakażeń.)

No więc chcę wydłubać dane postaci 92-K z tweetów publikowanych przez Ministerstwo Zdrowia. W tym celu za pomocą tweepy pobieram cały streamline (+3200 tweetów zaczynających się jeszcze w 2019 roku), ale dalej to już działam w Perlu, bo w Pythonie jakby mniej komfortowo się czuję. Po pierwsze zamieniam format json na csv:

use JSON;
use Data::Dumper;
use Time::Piece;
use open ":encoding(utf8)";
use open IN => ":encoding(utf8)", OUT => ":utf8";
binmode(STDOUT, ":utf8");

##  ID = tweeta ID; date = data; 
##  repid -- odpowiedź na tweeta o numerze ID
##  text -- tekst tweeta
print "id;date;repid;text\n";

while (<>) {
  chomp();
  $tweet =  $_;

  my $json = decode_json( $tweet );
  #print Dumper($json);
  $tid = $json->{"id"};
  $dat = $json->{"created_at"};
  ## Data jest w formacie rozwlekłym zamieniamy na YYYY-MM-DDTHH:MM:SS
  ## Fri Oct 04 14:48:25 +0000 2019
  $dat = Time::Piece->strptime($dat,
     "%a %b %d %H:%M:%S %z %Y")->strftime("%Y-%m-%dT%H:%M:%S");
  $rep = $json->{"in_reply_to_status_id"};
  $ttx = $json->{"full_text"}; $ttx =~ s/\n/ /g;
  ## Zamieniamy ; na , w tekście żeby użyć ; jako separatora
  $ttx =~ s/;/,/g; ####

  print "$tid;$dat;$rep;$ttx\n";

Komunikaty dłuższe niż limit Twittera są dzielone na kawałki, z których każdy jest odpowiedzią na poprzedni, np:

1298900644267544576;2020-08-27T08:30:20;1298900642522685440;53-M, 78-K i 84-K Kraków. Większość osób ...
1298900642522685440;2020-08-27T08:30:20;1298900640714948608;67-K Lublin (mieszkanka woj. podkarpackiego), 85-K Łańcut,...
1298900640714948608;2020-08-27T08:30:19;1298900639586680833;kujawsko-pomorskiego (24), świętokrzyskiego (18), opolskiego...
1298900639586680833;2020-08-27T08:30:19;;Mamy 887 nowych i potwierdzonych przypadków zakażenia #koronawirus z województw: ...

Czyli tweet 1298900639586680833 zaczyna, 1298900640714948608 jest odpowiedzią na 1298900639586680833, a 1298900642522685440 odpowiedzią na 1298900640714948608 itd. Tak to sobie wymyślili... W sumie chyba niepotrzebnie, ale w pierwszym kroku agreguję podzielone komunikaty w ten sposób, że wszystkie odpowiedzi są dołączane do pierwszego tweeta (tego z pustym polem in_reply_to_status_id):

## nextRef jest rekurencyjna zwraca numer-tweeta,
## który jest początkiem wątku
sub nextRef {
  my $i = shift;

  if ( $RR{"$i"} > 0 ) {  
    return ( nextRef( "$RR{$i}" ) );
  } else { return "$i" }
}

### ### ###
while (<>) { chomp();
   ($id, $d, $r, $t) = split /;/, $_;

   $TT{$id} = $t;
   $RR{$id} = $r;
   $DD{$id} = $d;
}

### ### ###
for $id ( sort keys %TT ) {  
   $lastId = nextRef("$id");
   $LL{"$id"} = $lastId;
   $LLIds{"$lastId"} = "$lastId"; 
}

### ### ###
for $id (sort keys %TT) {  
    ## print "### $DD{$id};$id;$LL{$id};$TT{$id}\n";  }
    $TTX{$LL{"$id"}} .= " " . $TT{"$id"};
    $DDX{$LL{"$id"}} .= ";" . $DD{"$id"};
}
### ### ###

for $i (sort keys %TTX) { 

    $dates = $DDX{"$i"};
    $dates =~ s/^;//; ## pierwszy ; jest nadmiarowy
    @tmpDat = split /;/, $dates;

    $dat_time_ = $tmpDat[0]; 
    ($dat_, $time_) = split /T/, $dat_time_;
    $ffN = $#tmpDat + 1;
    $collapsedTweet = $TTX{$i};
    print "$i;$dat_;$time_;$ffN;$collapsedTweet\n";
}

Zapuszczenie powyższego powoduje konsolidację wątków, tj. np. powyższe 4 tweety z 2020-08-27 połączyły się w jeden:

1298900639586680833;2020-08-27;08:30:19;4; Mamy 887 nowych i potwierdzonych przypadków
zakażenia #koronawirus z województw: małopolskiego (233), śląskiego (118), mazowieckiego (107),
[...] Liczba zakażonych koronawirusem: 64 689 /2 010 (wszystkie pozytywne przypadki/w tym osoby zmarłe).

Teraz wydłubujemy tylko tweety z frazą nowych i potwierdzonych albo nowe i potwierdzone:

## nowych i potwierdzonych albo nowe i potwierdzone
## (MZ_09.csv to CSV ze `skonsolidowanymi' wątkami)
cat MZ_09.csv | grep 'nowych i pot\|nowe i pot' > MZ_09_C19.csv
wc -l MZ_09_C19.csv
189 MZ_09_C19.csv

Wydłubanie fraz wiek-płeć ze skonsolidowanych wątków jest teraz proste:

  perl -e '
while (<>) {
 ($i, $d, $t, $c, $t) =  split /;/, $_;

  while ($t =~ m/([0-9]+-[MK])/g ) {
   ($w, $p) = split /\-/, $1;
   print "$d;$w;$p\n";
  }
}' MZ_09_C19.csv > C19D.csv
wc -l C19PL_down.csv
1738

Plik wykazał 1738 osób. Pierwsza komunikat jest z 16 kwietnia. Ostatni z 31. sierpnia. Pierwsze zarejestrowane zgony w PL odnotowano 12 marca (albo 13 nieważne). W okresie 12 marca --15 kwietnia zmarło 286 osób. Dodając 286 do 1738 wychodzi 2024. Wg MZ w okresie 12.03--31.08 zmarło 2039. Czyli manko wielkości 15 zgonów (około 0,5%). No trudno, nie chce mi się dociekać kto i kiedy pogubił tych 15...

Równie prostym skryptem zamieniam dane indywidualne na tygodniowe

#!/usr/bin/perl -w
use Date::Calc qw(Week_Number);

while (<>) {
  chomp();
  ##if ($_ =~ /age/) { next }  ##

  my ($d, $w, $p ) = split /;/, $_;
  my ($y_, $m_, $d_) = split /\-/, $d;
  my $week = Week_Number($y_, $m_, $d_);

  $DW{$week} += $w;
  $DN{$week}++;

  $DD{$d} = $week;
  $YY{$week} = 0;
  $YY_last_day{$week} = $d;

  ## wiek wg płci
  $PW{$p} += $w;
  $PN{$p}++;
}

for $d (keys %DD) {
    $YY{"$DD{$d}"}++; ## ile dni w tygodniu   
 }

print STDERR "Wg płci/wieku (ogółem)\n";

for $p (keys %PW) {
  $s = $PW{$p}/$PN{$p};
  printf STDERR "%s %.2f %i\n", $p, $s, $PN{$p};
  $total += $PN{$p};
}

print STDERR "Razem: $total\n";
print "week;deaths;age;days;date\n";

for $d (sort keys %DW) {
  if ($YY{$d} > 2 ) {## co najmniej 3 dni 
    $s = $DW{$d}/$DN{$d};
    printf "%s;%i;%.2f;%i;%s\n", $d, $DN{$d}, $s, $YY{$d}, $YY_last_day{$d};
  }
}

Co daje ostatecznie (week -- numer tygodnia w roku; meanage -- średni wiek zmarłych; deaths -- liczba zmarłych w tygodniu; days -- dni w tygodniu; date -- ostatni dzień tygodnia):

week;deaths;meanage;days;date
16;55;77.07;4;2020-04-19
17;172;75.09;7;2020-04-26
18;144;77.29;7;2020-05-03
19;123;76.46;7;2020-05-10
20;126;76.40;7;2020-05-17
21;71;76.37;7;2020-05-24
22;68;78.12;7;2020-05-31
23;93;75.73;7;2020-06-07
24;91;75.93;7;2020-06-14
25;109;77.24;7;2020-06-21
26;83;75.06;7;2020-06-28
27;77;74.09;7;2020-07-05
28;55;76.91;7;2020-07-12
29;54;77.33;7;2020-07-19
30;48;76.52;7;2020-07-26
31;60;74.88;7;2020-08-02
32;76;77.17;7;2020-08-09
33;71;73.11;7;2020-08-16
34;77;75.61;7;2020-08-23
35;79;74.33;7;2020-08-30

To już można na wykresie przedstawić:-)

library("dplyr")
library("ggplot2")
library("scales")
##
spanV <- 0.25
d <- read.csv("C19D_weekly.csv", sep = ';',  header=T, na.string="NA")

first <- first(d$date)
last <- last(d$date)
period <- sprintf ("%s--%s", first, last)
d$deaths.dailymean <- d$deaths/d$days

cases <- sum(d$deaths);
max.cases <- max(d$deaths)

note <- sprintf ("N: %i (source: twitter.com/MZ_GOV_PL)", cases)

pf <- ggplot(d, aes(x= as.Date(date), y=meanage)) +
 geom_bar(position="dodge", stat="identity", fill="steelblue") +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80)) +
 xlab(label="week") +
 ylab(label="mean age") +
 ggtitle(sprintf ("Mean age of COVID19 fatalities/Poland/%s", period), subtitle=note ) 

note <- sprintf ("N: %i (source: twitter.com/MZ_GOV_PL)", cases)
pg <- ggplot(d, aes(x= as.Date(date), y=deaths.dailymean)) +
 geom_bar(position="dodge", stat="identity", fill="steelblue") +
 scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
 scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40,45,50)) +
 xlab(label="week") +
 ylab(label="daily mean") +
 ggtitle(sprintf("Daily mean number of COVID19 fatalities/Poland/%s", period), subtitle=note )

ggsave(plot=pf, file="c19dMA_w.png")
ggsave(plot=pg, file="c19dN_w.png")

Wynik tutaj:

Średnia Wg płci/wieku/ogółem. Ogółem -- 75,9 lat (1738 zgonów) w tym: mężczyźni (914 zgonów) -- 73.9 lat oraz kobiety (824) -- 78.4 lat. Stan pandemii na 31.08 przypominam. Apogeum 2 fali... Średnia wieku w PL (2019 rok) 74,1 lat (M) oraz 81,8 lat (K).

Jak już pisałem danych nt COVID19 jest multum bo są traktowane jako treść promocyjna, przyciągająca klikających. Każda tuba medialna (gazeta/portal/telewizja) w szczególności publikuje dane nt.

Źródłem pierwotnym każdego wydają się być raporty narodowe (bo jak inaczej), ale ponieważ te raporty narodowe są składane w różny sposób, to ich połączenie w jedną bazę też różnie może wyglądać. Generalnie ci co takie bazy robią to albo przyznają się, że działają na hmmm niekonwencjonalnych źródłach (Twitter) albo nic nie piszą, skąd mają dane. Mają i już...

Wydaje się (chyba, że czegoś nie wiem), że ECDC, OWiD, CSSE oraz Worldometers (dalej WMs) są najpopularniejszymi źródłami danych nt COVID19 w przekroju międzynarodowym. (Nawiasem mówiąc: WHO nie publikuje danych -- publikuje raporty z danymi w formacie PDF. Wydobycie z nich danych jest nietrywialne i kosztowne, bo nie da się tego na 100% zautomatyzować. W rezultacie prawie nikt nie powołuje się na WHO jako źródło danych -- lekki szejm przyznajmy, bo niby ta organizacja jest od tego żeby m.in. zbierać i udostępniać informację n/t.) Taka drobna różnica na początek: ECDC, OWiD oraz CSSE to prawdziwe bazy: zarejestrowane z dzienną częstotliwością zgony, przypadki, testy i co tam jeszcze. OWiD kopiuje dane z ECDC, kiedyś kopiowało z WHO ale napisali że WHO zawierało liczne błędy i to ich skłoniło do korzystania z ECDC (0:2 dla WHO). WMs publikuje stan na, bazy jako takiej nie ma (przynajmniej publicznie albo nie potrafię jej odszukać na stronie). Można założyć że jak się ogląda stronę WMs z ,,notowaniami'' nt/ koronawirusa w dniu X o godzinie T to jest to stan na X:T. Nawiasem mówiąc tak jest wygodniej, ale jednocześnie komplikuje to sprawę w aspekcie: dzienna liczba przypadków chociażby z uwagi na różnice czasu (jak w PL kończy się dzień to na Fiji jest w połowie inny; inna sprawa, że wątpię żeby ktoś się tym przejmował). Niemniej WMs ma rubrykę "nowe przypadki", tyle że nie bardzo wiadomo co to znaczy...

No więc po tym przydługim wstępie do rzeczy: jak się mają dane z WMs względem ECDC? Jak wspomniałem, na stronie WMs nie ma bazy -- jest tabela z danymi ze stanem ,,na teraz''. ECDC z kolei publikuje bazę w postaci arkusza kalkulacyjnego. Ściągam dane codziennie. Ze strony WMs o 21:00 (koniec dnia, przynajmniej w PL) oraz o 23:00 ze strony ECDC. Dane te wyglądają jakoś tak (WMs, po konwersji HTML→CSV):

date;country;totalC;newC;totalD;newD;totalT
04040600;USA;277467;+306;7402;+10;830095

Stempel czasu jest ustalany w momencie pobrania danych. Na stronie WMs czas nie jest podany explicite (nie ma czegoś takiego jak np. dane aktualizowano o H:M). Czyli 04040600 to dane z 2020/04/04 z godziny 6:00.

Dane ECDC wyglądają jakoś tak:

date;id;country;newc;newd;totalc;totald
2020-04-04;US;United_States_of_America;32425;1104;277965;7157

NewC -- nowe przypadki (dzienne); NewD -- nowe zgodny; totalC -- przypadki łącznie; totalD -- zgony łącznie. Baza ECDC ma stempel czasu (dzień).

W przypadku PL wiem, że Ministerstwo Zdrowia (MinZ) publikuje dane generalnie o godzinie 10-coś-tam oraz o 17/18-coś-tam. (Czemu tak nie wiem). Patrząc na dane z WMs wiedzę, że o 21:00 publikują już dane aktualne na ten dzień, w tym sensie, że uwzględnią stan z ostatniego dziennego komunikatu MinZ (ale jakiego formalnie dnia te dane dotyczą, to już inna sprawa, bo ten dzień przecież się nie skończył :-)). Jeżeli chodzi o ECDC to dane pobrane w dniu X zawierają dane do dnia X-1, żeby było śmieszniej ECDC dla tego dnia przypisuje dane z komunikatu MinZ z dnia X-2. Czyli na przykładzie: arkusz pobrany o 23:00 dnia 24/04/2020 będzie miał ostatni wiersz datowany 23/04 ale dane w tym wierszu będą tymi które pojawiły się na stronie MinZ w dniu 22/04.

Uzbrojony o tę wiedzę dla wybranych 24 krajów wykreśliłem dane (z kwietnia) w wersji WMs oraz ECDC, w dwóch wariantach: z oryginalnymi stemplami czasowymi (górny wiersz) oraz ze stemplem skorygowanym przy założeniu że dane ECDC są 24H opóźnione (czyli dzień 23/04 tak naprawdę to dzień 22/04 itd). Te ,,skorygowane dane'' to dolny wiersz. Dla 90% krajów dane łącznie nakładają się czyli dane są identyczne (są wyjątki--ciekawe czemu). Dane dzienne to misz-masz, każda baza ma własną wersję, nie wiadomo która jest prawdziwa, tyle, że ECDC ma zawsze dane dzienne a WMs niekoniecznie (dla Japonii prawie zawsze ta kolumna była pusta)

Dane i komplet wykresów jest tutaj

Poniżej kilka wybranych krajów:

Marimekko to zestawiony do 100% wykres słupkowy, gdzie szerokość słupka jest proporcjonalna do jego udziału w liczebności. (https://predictivesolutions.pl/wykres-marimekko-czyli-analityczny-patchwork)

Nie wiem o co chodzi, ale patrząc na przykłady, to jest to wykres pokazujący strukturę dwóch zmiennych na raz. Coś jak skumulowany wykres słupkowy (stacked barchart), tyle że słupki mają zmienną szerokość, odpowiadającą udziałom/liczebnościom wartości jednej cechy. Dla każdego słupka z kolei poszczególne segmenty mają wysokości proporcjonalne do udziałów/liczebności wartości drugiej cechy (w tym słupku). Alternatywną nazwą jest wykres mozaikowy.

Ale jest też trochę inny wariant takich wykresów, dla przypadku kiedy dla każdej jednostki w populacji generalnej jest Wi = Ci/Ni. Jeżeli wysokość słupka jest proporcjonalna do wartości Wi, szerokość jest proporcjonalna do Ni, to pola są oczywiście w proporcji Ci. Przykładowo jeżeli populacją generalną są kraje świata, wartością cechy Ni jest liczba mieszkańców w kraju i, wartością cechy Ci wielkość emisja CO2 w kraju i, to Wi jest oczywiście emisją per capita.

Moim zdaniem użyteczne, bo pokazuje na raz dwa natężenia: łączne, w skali całej populacji oraz szczegółowe, w skali jednostki. Kontynuując przykład: ile emituje przeciętny mieszkaniec kraju i oraz jaki jest udział emisji kraju i w całości emisji.

W przypadku epidemii COVID19 podstawową zmienną jest liczba zarażonych/zmarłych w kraju i. Ale jeżeli chcemy porównać kraj i z krajem j to oczywiście należy uwzględnić liczbę mieszkańców w obu krajach. Czyli wysokości słupków powinny odpowiadać liczbie zarażonych/zmarłych na jednostkę (np. na 1mln) a szerokości liczbie mieszkańców:

library(ggplot2)
library(dplyr)

dat <- "2020/04/09"

d <- read.csv("indcs.csv", sep = ';',  header=T, na.string="NA");
## liczba ludności w milionach (szerokość):
d$popm <- d$pop / million

## Oblicz współczynniki na 1mln
d$casesr <- d$cases/d$popm
### Wysokość:
d$deathsr <- d$deaths/d$popm

## Ograniczamy liczbę krajów żeby zwiększyć czytelność wykresu
## Tylko kraje wykazujące zmarłych
d <- d %>% filter(deaths > 0) %>% as.data.frame
## Tylko kraje z min 2/1mln i populacji > 1mln
d9 <- d %>% filter(deathsr > 2 & popm > 3 & deaths > 49) %>% droplevels() %>%
 arrange (deathsr) %>% as.data.frame

d9$w <- cumsum(d9$popm)
d9$wm <- d9$w - d9$popm
d9$wt <- with(d9, wm + (w - wm)/2)

d8$w <- cumsum(d8$popm)
d8$wm <- d8$w - d8$popm
d8$wt <- with(d8, wm + (w - wm)/2)

## Dzielimy etykiety na dwie grupy
## (inaczej wiele etykiet zachodzi na siebie)
d9$iso3h <- d9$iso3
d9$iso3l <- d9$iso3
## Kraje o niskich  wartościach bez etykiet
d9$iso3h[ (d9$popm < 15 ) ] <- ""
## Kraje o wysokich  wartościach bez etykiet
d9$iso3l[ (d9$popm >= 15 ) ] <- ""

p9  <- ggplot(d9, aes(ymin = 0)) +
  ylab(label="mratio (deaths/1mln)") +
  xlab(label="population (mln)") +
  ggtitle(sprintf("COVID19 mortality (%s | mratio > 2 | population > 3mln )", dat), 
      subtitle="source: https://www.ecdc.europa.eu/en/covid-19-pandemic (twitter.com/tprzechlewski)") +
  geom_rect(aes(xmin = wm, xmax = w, ymax = deathsr, fill = iso3)) +
  geom_text(aes(x = wt, y = 0, label = iso3h), vjust=+0.5, hjust=+1.25, size=2.0, angle = 90) +
  geom_text(aes(x = wt, y = 0, label = iso3l), vjust=+0.5, hjust=-0.20, size=1.5, angle = 90) +
  theme(legend.position = "none") 
  ## ... podobnie dla zmiennej casesr

Wynik w postaci rysunków:

Powyższe jest dostępne tutaj

Wymyśliłem sobie wykreślić wykresy pokazujące zależność pomiędzy liczbą zarażonych/zmarłych, a wybranymi wskaźnikami: GDP (zamożność), oczekiwana długość życia (poziom służby zdrowia) oraz śmiertelność dzieci (zamożność/poziom rozwoju). Większość danych pobrałem z portalu OWiD:

• GDP (https://ourworldindata.org/economic-growth maddison-data-gdp-per-capita-in-2011us.csv)
• LE (https://ourworldindata.org/life-expectancy life-expectancy.csv)
• CM (https://ourworldindata.org/child-mortality child-mortality-igme.csv).

Dane dotyczące liczby ludności są z bazy Banku Światowego (https://data.worldbank.org/indicator/sp.pop.totl.) Dane dotyczące zarażonych/zmarłych ze strony ECDC (www.ecdc.europa.eu/en/covid-19/data-collection) Scaliłem wszystko do kupy skryptem Perlowym. NB pojawił się tzw. small problem, bo bazy OWiD oraz WorldBank używają ISO-kodów 3-literowych krajów, a ECDC dwuliterowych (PL vs POL). Trzeba było znaleźć wspólny mianownik. Szczęśliwie Perl ma gotowy moduł. Oprócz tego ECDC stosuje EU-standard w przypadku Grecji (EL zamiast GR) oraz Wielkiej Brytanii (UK/GB).

#!/usr/bin/perl
use Locale::Codes::Country;
...
while (<COVID>) { chomp();
($date, $iso2, $country, $newc, $newd, $totalc, $totald) = split /;/, $_;

   $iso3 = uc ( country_code2code($iso2, 'alpha-2', 'alpha-3'));
}

Rezultat jest zapisywany do pliku o następującej zawartości:

iso3;country;lex2019;gdp2016;cm2017;pop2018;cases;deaths
ABW;Aruba;76.29;NA;NA;105845;77;0
AFG;Afghanistan;64.83;1929;6.79;37172386;423;14
...

Wierszy jest 204, przy czym niektórym krajom brakuje wartości. Ponieważ dotyczy to krajów egzotycznych, zwykle małych, to pominę je (nie teraz później, na etapie przetwarzania R-em). Takich wybrakowanych krajów jest 49 (można grep-em sprawdzić). Jedynym większym w tej grupie jest Syria, ale ona odpada z innych powodów.

Do wizualizacji wykorzystam wykres punktowy (dot-plot) oraz wykresy rozrzutu (dot-plot).

library("dplyr")
library("ggplot2")
library("ggpubr")
## 
options(scipen=1000000)
## https://www.r-bloggers.com/the-notin-operator/
`%notin%` <- Negate(`%in%`)
##
today <- Sys.Date()
tt<- format(today, "%d/%m/%Y")
million <- 1000000
## Lista krajów Europejskich + Izrael
## (pomijamy kraje-liliputy)
ee <- c(
'BEL', 'GRC', 'LTU', 'PRT', 'BGR', 'ESP', 'LUX', 'ROU', 'CZE', 'FRA', 'HUN',
'SVN', 'DNK', 'HRV', 'MLT', 'SVK', 'DEU', 'ITA', 'NLD', 'FIN', 'EST', 'CYP',
'AUT', 'SWE', 'IRL', 'LVA', 'POL', 'ISL', 'NOR', 'LIE', 'CHE', 'MNE', 'MKD',
'ALB', 'SRB', 'TUR', 'BIH', 'BLR', 'MDA', 'UKR', 'ISR', 'RUS', 'GBR' );
ee.ee <- c('POL')

d <- read.csv("indcs.csv", sep = ';',  header=T, na.string="NA");
## Liczba krajów
N1 <- nrow(d)
## liczba ludności w milionach
d$popm <- d$pop / million

## Oblicz współczynniki na 1mln 
d$casesr <- d$cases/d$popm 
d$deathsr <- d$deaths/d$popm 

## Tylko kraje wykazujące zmarłych
d <- d %>% filter(deaths > 0) %>% as.data.frame
## Liczba krajów wykazujących zmarłych
N1d <- nrow(d)

## Tylko kraje z kompletem wskaźników (pomijamy te z brakami)
d <- d[complete.cases(d), ]

nrow(d)

##  UWAGA: pomijamy kraje o wsp. <= 2
##  droplevels() usuwa `nieużywane' czynniki
##  mutate zmienia kolejność na kolejność wg dearhsr:
d9 <- d %>% filter(deathsr > 2 ) %>% droplevels() %>% 
 mutate (iso3 = reorder(iso3, deathsr)) %>% as.data.frame

N1d2 <- nrow(d9)
M1d2 <- median(d9$deathsr, na.rm=T)

## https://stackoverflow.com/questions/11093248/geom-vline-with-character-xintercept
## Wykres punktowy
rys99 <- ggplot(d9, aes(x =iso3, y = deathsr )) +
  geom_point(size=1, colour = 'steelblue', alpha=.5) +
  xlab(label="Country") +
  ylab(label="Deaths/1mln") +
  ggtitle(sprintf("COVID19 mortality in deaths/mln (as of %s)", tt), 
   subtitle=sprintf("Countries with ratio > 0: %i | Countries with ratio > 2.0: N=%i (median %.1f)", 
       N1d, N1d2, M1d2)) +
  theme(axis.text = element_text(size = 4)) +
  ##theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(breaks=c(0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300,320,340,360)) +
  geom_hline(yintercept=M1d2, linetype="solid", color = "steelblue") +
  geom_vline(aes(xintercept = which(levels(iso3) == 'POL')), size=1, color="#8A0303", alpha=.25) +
  geom_vline(aes(xintercept = which(levels(iso3) == 'DEU')), size=1, color="#8A0303", alpha=.25) +
  geom_vline(aes(xintercept = which(levels(iso3) == 'SWE')), size=1, color="#8A0303", alpha=.25) +
  coord_flip()

Uwaga: oś OX to skala porządkowa. Czynniki powinny być uporządkowane wg. wartości zmiennej z osi OY, czyli według deathsr. Do tego służy funkcja mutate (iso3 = reorder(iso3, deathsr). Można też uporządkować je ,,w locie'' aes(x =iso3, y = deathsr ), ale wtedy niepoprawnie będą kreślone (za pomocą geom_vline) linie pionowe. Linie pionowe służą do wyróżnienia pewnych krajów. Linia pozioma to linia mediany.

sources <- sprintf ("As of %s\n(Sources: %s %s)", tt,
  "https://www.ecdc.europa.eu/en/covid-19-pandemic",
  "https://ourworldindata.org/")

## Żeby etykiety nie zachodziły na siebie tylko dla wybranych krajów
## Add empty factor level! Istotne inaczej będzie błąd
# https://rpubs.com/Mentors_Ubiqum/Add_Levels
d$iso3 <- factor(d$iso3, levels = c(levels(d$iso3), ""))

d$iso3xgdp <- d$iso3
d$iso3xlex <- d$iso3
d$iso3xcm <- d$iso3

## Kraje o niskich  wartościach bez etykiet
## Bez etykiet jeżeli GDP<=45 tys oraz wskaźnik < 50:
d$iso3xgdp[ (d$gdp2016 < 45000) & (d$deathsr < 50 ) ] <- ""
## Inne podobnie
d$iso3xlex[ ( (d$lex2019 < 80) | (d$deathsr < 50 ) ) ] <- ""
d$iso3xcm[ ((d$cm2017 > 1.2) | (d$deathsr < 50 ) ) ] <- ""

## GDP vs współczynnik zgonów/1mln
rys1 <- ggplot(d, aes(x=gdp2016, y=deathsr)) + 
  geom_point() +
  geom_text(data=d, aes(label=sprintf("%s", iso3xgdp), x=gdp2016, y= deathsr), vjust=-0.9, size=2 ) +
  xlab("GDP (USD, Constant prices)") + 
  ylab("deaths/1mln") + 
  geom_smooth(method="loess", se=F, size=2) +
  ggtitle("GDP2016CP vs COVID19 mortality", subtitle=sources)

## Life ex vs współczynnik zgonów/1mln
rys2 <- ggplot(d, aes(x=lex2019, y=deathsr)) + 
  geom_point() +
  geom_text(data=d, aes(label=sprintf("%s", iso3xlex), x=lex2019, y= deathsr), vjust=-0.9, size=2 ) +
  xlab("Life expentancy") + 
  ylab("deaths/1mln") + 
  geom_smooth(method="loess", se=F, size=2) +
  ggtitle("Life expentancy vs COVID19 mortality", subtitle=sources)

## Child mortality vs współczynnik zgonów/1mln
rys3 <- ggplot(d, aes(x=cm2017, y=deathsr)) + 
  geom_point() +
  geom_text(data=d, aes(label=sprintf("%s", iso3xcm), x=cm2017, y= deathsr), vjust=-0.9, size=2 ) +
  xlab("Child mortality %") +
  ylab("deaths/1mln") + 
  geom_smooth(method="loess", se=F, size=2) +
  ggtitle("Child mortality vs COVID19 mortality", subtitle=sources)

## GDP vs Child mortality 
rys0 <- ggplot(d, aes(x=gdp2016, y=cm2017)) + 
  geom_point() +
  xlab("GDP (USD, Constant prices)") + 
  ylab("Child mortality") +
  geom_smooth(method="loess", se=F, size=2) +
  ggtitle("GDP2016CP vs Child mortality", subtitle=sources)

Jeszcze raz dla krajów Europejskich:

## Tylko kraje Europejskie:
d <- d %>% filter (iso3 %in% ee) %>% as.data.frame

d$iso3xgdp <- d$gdp2016
d$iso3xlex <- d$lex2019
d$iso3xcm <- d$cm2017
d$iso3xgdp[ d$iso3 %notin% ee.ee ] <- NA
d$iso3xlex[ d$iso3 %notin% ee.ee ] <- NA
d$iso3xcm[ d$iso3 %notin% ee.ee ] <- NA

rys1ec <- ggplot(d, aes(x=gdp2016, y=casesr)) + 
  geom_point() +
  geom_text(data=d, aes(label=sprintf("%s", iso3), x=gdp2016, y= casesr), vjust=-0.9, size=2 ) +
  geom_point(data=d, aes(x=iso3xgdp, y= casesr), size=2, color="red" ) +
  xlab("GDP (USD, Constant prices") + 
  ylab("cases/1mln") +
  geom_smooth(method="loess", se=F, size=2) +
  ggtitle("GDP2016CP vs COVID19 cases (Europe)", subtitle=sources)
  ... itd...

Wynik w postaci rysunków:

Powyższe jest dostępne tutaj

Google has launched a new website that uses anonymous location data collected from users of Google products and services to show the level of social distancing taking place in various locations. The COVID-19 Community Mobility Reports web site will show population data trends of six categories: Retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residential. The data will track changes over the course of several weeks, and as recent as 48-to-72 hours prior, and will initially cover 131 countries as well as individual counties within certain states. (cf. www.google.com/covid19/mobility/.)

The raports contains charts and comments in the form: NN% compared to baseline (in six above mentioned categories) where NN is a number. It is assumed the number is a percent change at the last date depicted (which accidentaly is a part of a filename). So for example a filename 2020-03-29_PL_Mobility_Report_en.pdf contains a sentence `Retail & recreation -78% compared to baseline` which (probably) means that (somehow) registered traffic at R&R facilities was 22% of the baseline. Anyway those six numbers was extracted for OECD countries (and some other countries) and converted to CSV file.

The conversion was as follows: first PDF files was downloaded with simple Perl script:

#!/usr/bin/perl
# https://www.google.com/covid19/mobility/
use LWP::UserAgent;
use POSIX 'strftime';

my $sleepTime = 11;

%OECD = ('Australia' => 'AU', 'New Zealand' => 'NZ',
'Austria' => 'AT', 'Norway' => 'NO', 'Belgium' => 'BE',
'Poland' => 'PL', 'Canada' => 'CA', 'Portugal' => 'PT',
'Chile' => 'CL', 'Slovak Republic' => 'SK',
## etc ...
);

@oecd = values %OECD;

my $ua = LWP::UserAgent->new(agent => 'Mozilla/5.0', cookie_jar =>{});
my $date = "2020-03-29";

foreach $c (sort @oecd) {
   $PP="https://www.gstatic.com/covid19/mobility/${date}_${c}_Mobility_Report_en.pdf";

   my $req = HTTP::Request->new(GET => $PP);
   my $res = $ua->request($req, "${date}_${c}_Mobility_Report_en.pdf");

   if ($res->is_success) { print $res->as_string; }
   else { print "Failed: ", $res->status_line, "\n"; }
}

Next PDF files was converted to .txt with pdftotext. The relevant fragments of .txt files looks like:

  Retail & recreation
+80%

-78%
compared to baseline

So it looks easy to extract the relevant numbers: scan line-by-line looking for a line with appropriate content (Retail & recreation for example). If found start searching for 'compared to baseline'. If found retrieve previous line:

#!/usr/bin/perl
$file = $ARGV[0];

while (<>) {   chomp();
  if (/Retail \& recreation/ ) { $rr = scan2base(); }
  if (/Grocery \& pharmacy/ ) { $gp = scan2base(); }
  if (/Parks/ ) { $parks = scan2base(); }
  if (/Transit stations/ ) { $ts = scan2base(); }
  if (/Workplaces/ ) { $wps = scan2base(); }
  if (/Residential/ ) { $res = scan2base();
     print "$file;$rr;$gp;$parks;$ts;$wps;$res\n";
     last;  }
}

sub scan2base {
  while (<>) {
   chomp();
   if (/compared to baseline/) {  return ($prevline); }
   $prevline = $_;
  }
}

Extracted data can be found here.

Danych nt COVID19 jest multum bo są traktowane jako treść promocyjna, przyciągająca klikających. Każda tuba medialna (gazeta/portal/telewizja) w szczególności publikuje dane nt

Niestety z faktu, że danych nt COVID19 jest multum niewiele wynika, bo wszystkie są do dupy, w sensie że są wątpliwej jakości, tj. zwykle sposób w jaki są gromadzone nie jest opisany. Nie wiadomo kto jest klasyfikowany jako zarażony COVID19, nie wiadomo kto jest klasyfikowany jako zmarły w wyniku zarażenia COVID19. Można się domyślać że klasyfikowany jako zarażony COVID19 jest ten komu wykonany słynny test (w większości wypadków, podobno nie zawsze); zmarły w wyniku zarażenia COVID19 jest ten, któremu lekarz wypisał świadectwo zgonu ze stosownym wpisem.

Powyższe skutkuje: niemożnością oceny prawdziwej skali zjawiska (stąd teorie że rząd fałszuje) oraz niemożnością dokonania wiarygodnych porównań międzynarodowych.

Jeżeli chodzi o Polskę, to nikt nie prowadzi publicznego rejestru. Strona GIS to w ogóle kuriozalnie wygląda. Są komunikaty, jak ktoś ma czas to może jest sobie z nich dane wydłubywać i agregować. Na poziomie międzynarodowym są 2 źródła agregacji pierwotnej nazwijmy to: WHO oraz ECDC. Te dwa źródła agregują dane nadsyłane przez ciała krajowe, wg jakiejś niezdefiniowanej (przypuszczalnie ad hoc ustalanej) procedury. Inni korzystają z danych WHO/ECDC pośrednio lub bezpośrednio ewentualnie uzupełniając/modyfikując je w bliżej niezdefiniowany sposób. No i są jeszcze źródła specyficzne takie jak Google Community Mobility Reports.

WHO Situation Reports. To nie jest baza danych, ale pliki PDF zawierające raporty w tym dane. Pozyskanie z nich danych wymaga nietrywialnej konwersji. www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports . Dane z raportów dostępne są m.in. na stronie Wikipedii: en.wikipedia.org/wiki/2019%E2%80%9320_coronavirus_pandemic_cases/WHO_situation_reports oraz en.wikipedia.org/wiki/Talk:2019%E2%80%9320_coronavirus_pandemic_cases/WHO_situation_reports

ECDC.europa.eu Dane udostępniane w postacji codziennie aktualizowanego arkusza kalkulacyjnego. www.ecdc.europa.eu/en/covid-19/data-collection [Since the beginning of the coronavirus pandemic, ECDC's Epidemic Intelligence team has been collecting the number of COVID-19 cases and deaths, based on reports from health authorities worldwide.]

John Hopkins Univ/CSSE github.com/CSSEGISandData/COVID-19 [To identify new cases, we monitor various twitter feeds, online news services, and direct communication sent through the dashboard. Before manually updating the dashboard, we confirm the case numbers using regional and local health departments, namely the China CDC (CCDC), Hong Kong Department of Health, Macau Government, Taiwan CDC, European CDC (ECDC), the World Health Organization (WHO), as well as city and state level health authorities.]

Worldometers https://worldometers.info/coronavirus/ [nie wiadomo jak zbierane, przypuszczalnie kopiowane z WHO/ECDC; Worldometers, to -- wydaje się -- inicjatywa PR-owa firmy produkującej oprogramowanie]

OWiD czyli Our World in Data wykorzystuje bazę ECDC. ourworldindata.org/coronavirus-source-data [na podstawie ECDC]

Reasumując: jak ktoś potrzebuje gotowego zbioru danych, to ma do wyboru ECDC/OWiD/CSSE. Wszystkie są wątpliwe, ale lepszych nie ma a ci przynajmniej podają (ogólnikowo to fakt) jak te dane zbierają. Jak ktoś używa worldometers to pytanie czemu to robi... Jak posługuje się jeszcze innymi bardziej egoztycznymi danymi to szkoda tracić czasu na jego analizy (ew. sprawdzić czy nie są to dane ECDC/OWiD/CSSE tylko pod inną marką sprzedawane).

W Polsce nie ma oficjalnego rejestru. Przynajmniej ja nic nie wiem na temat. To tak nawiasem mówiąc szejm. Że żaden urząd, uniwersytet czy instytut nie udostępnia oficjalnych/wiarygodnych/kompletnych/łatwo dostępnych danych (w Niemczech na przykład robi to słynny RKI; a we Francji nie mniej słynny pasteur.fr). W PL zaś każdy się stara i coś tam udostępnia, z naciskiem na coś... Znalazłem rejestr nieopisany (w sensie jak/skąd są nim gromadzone dane) prowadzony przez dziennik z grupy PolskaPress. dziennikzachodni.carto.com/tables/zachorowania_na_koronawirusa_w_polsce_marzec/public

Google Community Mobility Reports To nie jest baza danych, ale zbiór raportów w formacie PDF. www.google.com/covid19/mobility/. [Google has launched a new website that uses anonymous location data collected from users of Google products and services to show the level of social distancing taking place in various locations. The COVID-19 Community Mobility Reports web site will show population data trends of six categories: Retail and recreation, grocery and pharmacy, parks, transit stations, workplaces, and residential. The data will track changes over the course of several weeks, and as recent as 48-to-72 hours prior, and will initially cover 131 countries as well as individual counties within certain states.] Ciekawostka raczej, bo w szczególności, nie do końca wiadomo co te procenty Gógla oznaczają, np. -60% względem baseline. Anie nie wiadomo co to jest ten baseline (średnia?) ani jak liczony jest ruch...

Nie mniej wydłubałem te procenty z raportów dla krajów OECD i zamieniłem na plik w formacie CSV. Jest on do pobrania tutaj.

Dane dotyczące USA. Oczywiście są częścią WHO/ECDC/CSSE. Ale są także bardziej szczegółowe:

CDC [The provisional counts for coronavirus disease (COVID-19) deaths are based on a current flow of mortality data in the National Vital Statistics System.] https://www.cdc.gov/nchs/nvss/vsrr/COVID19/index.htm

NewYork Times [The data is the product of dozens of journalists working across several time zones to monitor news conferences, analyze data releases and seek clarification from public officials on how they categorize cases.] https://github.com/nytimes/covid-19-data oraz https://www.nytimes.com/interactive/2020/us/coronavirus-us-cases.html

No i jeszcze są pewnie jakieś chińskie dane, ale to trzeba znać chiński.

Przez przypadek ciekawe odkrycie. W niedzielę zakończyły się wybory w Bawarii, wyłącznie w trybie korespondencyjnym z uwagi na epidemię #COVID19. Jednocześnie w PL trwa wałkowanie tematu pn przesunąć wybory prezydenckie. Niewątpliwie przykład niemiecki to kłopot, że tak powiem narracyjny, dla tych co chcą przesunięcia.

No więc naiwnie wpisałem w google: bavaria+second+round+elections+postal a w rezultacie dostałem głównie strony o apelu kandydatki Kidawy o przesunięcie wyborów w tym kuriozalna relacja Reutersa -- kiedyś szanowanej agencja informacyjnej, teraz kandyjskiej dezinformacyjnej prop-tuby. W dziale "Zdrowie" -- a jakże -- donosi ona o apelu o bojkot p. Kidawy kończąc ten "zdrowotny raport" raportem pana J. Flisa z Krakowa (tak to ten sam, udający naukowca, telewizyjny-profesor Flis), ale na temat wyborów w Bawarii, które skutkowały wg. p. Flisa 2 tys ofiar (metody wyliczeń, którą posłużył się "profesor" nie podano). Zaistne niezwykle relewantny dokument do mojego zapytania.

Financial Times zamieścił wykres wględnego tempa wzrostu (rate of growth) czyli procentu liczonego jako liczba-nowych / liczba-ogółem-z-okresu-poprzedniego x 100%. Na wykresie wględnego tempa wzrostu zachorowań na COVID19 wszystkim spada: Every day the Covid-19 virus is infecting an increasing number of people, but the rate of growth in cases in some of the worst-hit countries is starting to slow. Powyższe Czerscy przetłumaczyli jako m.in. trend dotyczy niemal wszystkich krajów rozwiniętych. [he, he... Rozwiniętych pod względem liczby chorych, pewnie chcieli uściślić, ale się nie zmieściło]

Spróbowałem narysować taki wykres samodzielnie:

library("dplyr")
library("ggplot2")
library("ggpubr")
##
surl <- "https://www.ecdc.europa.eu/en/publications-data/download-todays-data-geographic-distribution-covid-19-cases-worldwide"
today <- Sys.Date()
tt<- format(today, "%d/%m/%Y")

#d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", stringsAsFactors=FALSE);
d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", 
   colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));

d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)

Zwykłe read_csv skutkowało tym, że newc/newd nie były liczbami całkowitymi, tylko czynnikami. Z kolei dodanie colClasses kończyło się błędem. W końcu stanęło na tym, że czytam dane w kolumnach newc/newd zadeklarowanych jako napisy a potem konwertuję na liczby. Czy to jest prawidłowa strategia to ja nie wiem...

Kolejny problem: kolumny newc/newd zawierają NA, wykorzystywana później funkcja cumsum z pakietu dplyr, obliczająca szereg kumulowany nie działa poprawnie jeżeli szereg zawiera NA. Zamieniam od razu NA na zero. Alternatywnie można korzystać z funkcji replace_na (pakiet dplyr):

# change NA to 0
d[is.na(d)] = 0

# Alternatywnie replace_na
#d %>% replace_na(list(newc = 0, newd=0)) %>%
#  mutate( cc = cumsum(newc), dd=cumsum(newd))

Ograniczam się tylko do danych dla wybranych krajów, nie starszych niż 16 luty 2020:

d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame
str(d)

last.obs <- last(d$date)
c1 <- c('IT', 'DE', 'ES', 'UK', 'FR')
d1 <- d %>% filter (id %in% c1) %>% as.data.frame

str(d1)

Obliczam wartości skumulowane (d zawiera już skumulowane wartości, ale obliczone Perlem tak nawiasem mówiąc):

t1 <- d1 %>% group_by(id) %>%  summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T))

t1c <- d %>% group_by(id) %>%  mutate(cum_cc = cumsum(newc), cum_dd = cumsum(newd)) %>% 
  filter (id %in% c1) %>%
  filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame

  str(t1c)

Wykres wartości skumulowanych:

pc1c <- ggplot(t1c, aes(x= as.Date(date, format="%Y-%m-%d"), y=cum_cc)) + 
  geom_line(aes(group = id, color = id), size=.8) +
  xlab(label="") +
  theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
  ggtitle(sprintf("COVID19: total confirmed cases (%s)", last.obs), subtitle=sprintf("%s", surl)) 

ggsave(plot=pc1c, "Covid19_1c.png", width=15)

Kolumny cum_lcc/cum_ldd zawierają wartości z kolumny cum_cc/cum_dd ale opóźnione o jeden okres (funkcja lag):

## 
t1c <- t1c %>% group_by(id) %>% mutate(cum_lcc = lag(cum_cc)) %>% as.data.frame
t1c <- t1c %>% group_by(id) %>% mutate(cum_ldd = lag(cum_dd)) %>% as.data.frame

t1c$gr_cc <- t1c$newc / (t1c$cum_lcc + 0.01) * 100
t1c$gr_dd <- t1c$newd / (t1c$cum_ldd + 0.01) * 100

## Początkowo wartości mogą być ogromne zatem
## zamień na NA jeżeli gr_cc/dd > 90
t1c$gr_cc[ (t1c$gr_cc > 90) ] <- NA
t1c$gr_dd[ (t1c$gr_dd > 90) ] <- NA

Wykres tempa wzrostu:

pc1c_gr <- ggplot(t1c, aes(x= as.Date(date, format="%Y-%m-%d"), y=gr_cc,  colour = id, group=id )) + 
  ##geom_line(aes(group = id, color = id), size=.8) +
  geom_smooth(method = "loess", se=FALSE) +
  xlab(label="") +
  theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
  ggtitle(sprintf("COVID19: confirmed cases growth rate (smoothed)"), 
      subtitle=sprintf("%s", surl)) 

ggsave(plot=pc1c_gr, "Covid19_1g.png", width=15)

To samo co wyżej tylko dla PL/CZ/SK/HU:

c2 <- c('PL', 'CZ', 'SK', 'HU')

t2c <- d %>% group_by(id) %>%  mutate(cum_cc = cumsum(newc), cum_dd = cumsum(newd)) %>% 
  filter (id %in% c2) %>%
  filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame

##str(t2c)
t2c.PL <- t2c %>% filter (id == "PL") %>% as.data.frame
t2c.PL
head(t2c.PL, n=200)

pc2c <- ggplot(t2c, aes(x= as.Date(date, format="%Y-%m-%d"), y=cum_cc)) + 
  geom_line(aes(group = id, color = id), size=.8) +
  xlab(label="") +
  theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
  ggtitle(sprintf("COVID19: total confirmed cases (%s)", last.obs), subtitle=sprintf("Total: %s\n%s", lab1c, surl)) 

ggsave(plot=pc2c, "Covid19_2c.png", width=15)

t2c <- t2c %>% group_by(id) %>% mutate(cum_lcc = lag(cum_cc)) %>% as.data.frame
t2c <- t2c %>% group_by(id) %>% mutate(cum_ldd = lag(cum_dd)) %>% as.data.frame

t2c$gr_cc <- t2c$newc / (t2c$cum_lcc + 0.01) * 100
t2c$gr_dd <- t2c$newd / (t2c$cum_ldd + 0.01) * 100

## zamień na NA jeżeli gr_cc/dd > 90
t2c$gr_cc[ (t2c$gr_cc > 90) ] <- NA
t2c$gr_dd[ (t2c$gr_dd > 90) ] <- NA

t2c.PL <- t2c %>% filter (id == "PL") %>% as.data.frame
t2c.PL

pc2c_gr <- ggplot(t2c, aes(x= as.Date(date, format="%Y-%m-%d"), y=gr_cc,  colour = id, group=id )) + 
  ##geom_line(aes(group = id, color = id), size=.8) +
  geom_smooth(method = "loess", se=FALSE) +
  xlab(label="") +
  theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
  ggtitle(sprintf("COVID19: confirmed cases growth rate (smoothed)"), 
      subtitle=sprintf("%s", surl)) 

ggsave(plot=pc2c_gr, "Covid19_2g.png", width=15)

Koniec

Dane pierwotne: Center for Systems Science and Engineering (CSSE/Johns Hopkins University) https://github.com/CSSEGISandData/COVID-19 (także słynna wizualizacja: https://gisanddata.maps.arcgis.com/apps/opsdashboard/index.html#/bda7594740fd40299423467b48e9ecf6.) European Centre for Disease Prevention and Control https://ecdc.europa.eu/en/geographical-distribution-2019-ncov-cases.

Dane agregowane od innych: Our World in Data/Coronavirus Source Data/WHO Situation Reports https://ourworldindata.org/coronavirus-source-data . Są też tzw. dane w czasie rzeczywistym: https://worldometers.info/coronavirus/, ale ich wiarygodność jest podejrzana, bo w przeciwieństwie do tych wyżej opisanych nie wiadomo jak są zbierane i/lub skąd agregowane (więc nie ma klikalnego linku).

Na stronie https://ourworldindata.org/coronavirus-source-data są dane nt liczby przypadków/zgonów z powodu zarażenia wirusem covid19, których źródłem są 'Raporty Sytuacyjne WHO' (https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports/). Raporty są w formacie PDF więc bezpośrednio nie można korzystać z publikowanych tam danych. No ale uprzejmi ludzie z ourworldindata.org już te raporty zamienili na csv i są one gotowe do pobrania:

wget -N https://covid.ourworldindata.org/data/full_data.csv

Za pomocą prostych skryptów Perla modyfikuję plik full_data.csv tak, żeby poszczególne kolumny zawierały: date;id;country;newc;newd;totalc;totald (data, ISO-kod kraju, nazwa-kraju, nowe-przypadki, nowe-zgony, wszystkie-przypadki, wszystkie-zgony)

Dla wybranych krajów rysują wykresy liniowe (wykorzystując R):

library("dplyr")
library("ggplot2")
library("ggpubr")
##
today <- Sys.Date()
tt<- format(today, "%d/%m/%Y")

d <- read.csv("covid19.csv", sep = ';',  header=T, na.string="NA");

d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame

##
c1 <- c('ITw', 'DEw', 'ESw', 'UKw', 'FRw', 'DKw', 'SEw')
# date;id;country;newc;newd;totalc;totald
d1 <- d %>% filter (id %in% c1) %>% as.data.frame
t1 <- d1 %>% group_by(id) %>%  summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T))

lab1c <- toString(paste (sep=" = ", t1$id, t1$cc))
lab1d <- toString(paste (sep=" = ", t1$id, t1$dd))

str(d1)

pc1 <- ggplot(d1, aes(x= as.Date(date, format="%Y-%m-%d"), y=newc)) + geom_line(aes(group = id, color = id), size=.8) +
 xlab(label="") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle(sprintf("COVID19: new confirmed cases (%s)", tt), subtitle=sprintf("Total: %s\n%s", lab1c, surl))

pd1 <- ggplot(d1, aes(x= as.Date(date, format="%Y-%m-%d"), y=newd)) + geom_line(aes(group = id, color = id), size=.8) +
 xlab(label="") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle(sprintf ("COVID19: deaths (%s)", tt), subtitle=sprintf("Total: %s\n%s", lab1d, surl))

c2 <- c('PLw', 'CZw', 'SKw', 'HUw', 'ROw', 'BGw', 'ELw')
d2 <- d %>% filter (id %in% c2) %>% as.data.frame
t2 <- d2 %>% group_by(id) %>%  summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T))

str(d2)

lab2c <- toString(paste (sep=" = ", t2$id, t2$cc))
lab2d <- toString(paste (sep=" = ", t2$id, t2$dd))

pc2 <- ggplot(d2, aes(x= as.Date(date, format="%Y-%m-%d"), y=newc)) + geom_line(aes(group = id, color = id), size=.8) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 xlab(label="") +
 ggtitle(sprintf("COVID19: new confirmed cases (%s)", tt), subtitle=sprintf("Total: %s\n%s", lab2c, surl))

pd2 <- ggplot(d2, aes(x= as.Date(date, format="%Y-%m-%d"), y=newd)) + geom_line(aes(group = id, color = id), size=.8) +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 xlab(label="") +
 scale_y_continuous(breaks=c(1,2,3,4,5,6,7,8,9)) +
 ggtitle(sprintf ("COVID19: deaths (%s)", tt), subtitle=sprintf("Total: %s\n%s", lab2d, surl))

p1 <- ggarrange(pc1,pd1, ncol=2, nrow=1)
p2 <- ggarrange(pc2,pd2, ncol=2, nrow=1)
ggsave(plot=p1, "Covid19_1w.png", width=15)
ggsave(plot=p2, "Covid19_2w.png", width=15)

Zatem: Liczba przypadków/zgonów z powodu zarażenia wirusem covid19 na podstawie danych ourworldindata.org/WHO:

UE, która jest powszechnie krytykowana, że nic nie robi w sprawie, okazuje się że coś tobi -- też udostępnia jakieś dane w temacie (https://www.ecdc.europa.eu/en/publications-data/download-todays-data-geographic-distribution-covid-19-cases-worldwide). Wprawdzie nie jest określone skąd te dane pochodzą, ale sądząc po ich zawartości źródło jest to samo (WHO).

wget -N https://www.ecdc.europa.eu/sites/default/files/documents/COVID-19-geographic-disbtribution-worldwide-2020-03-15.xls\
  -O covid19.csv

Ponieważ dane, że tak powiem, unijne są w formacie xls zamieniam je na csv, wykorzystując do tego LibreOffice:

## zamień wszystkie pliki z bieżącego katalogu na csv ze ; (59) jako znakiem separacji:
soffice --convert-to csv:"Text - txt - csv (StarCalc)":59,,0,1,1 --outdir . *.xls

Rysuję wykresy liniowe zmodyfikowanym z dokładnością do pliku z danymi R-skryptem. Wyniki są prawie takie same. Może bym nawet nie zwrócił uwagi, że się różnią gdyby nie podejrzane załamanie liczby przypadków dla Włoch dla 15.03.2020 (z 2,5 tys na 90).

Drążąc temat wyrysowałem wykresy dla wybranych czterech krajów w dwóch wariantach danych (dane z ourworldindata.org oznaczone literką w). W szczególności i niestety Włochy 15/3/2020 odnotowały ponad 3497 nowych przypadków a nie 90 jak podano w bazie Unijnej. Są też mniejsze różnice w innych miejscach:

Wszystko to robione jest automatem co pobiera/zamienia/rysuje/wstawia na githuba (https://github.com/hrpunio/Nafisa/tree/master/Covid19) oraz wysyła na twittera (https://twitter.com/tprzechlewski). Automat działa na RaspberryPi zresztą...