Że czwarta fala szaleje to nadmiarowe zgony w Europie sobie liczę, w której to kategorii Polska cytując poetę trzecie miejsce w świecie (Kazik Staszewski/Amnezja)

Dane pobieram za pomocą get_eurostat (pakiet eurostat):

library("eurostat")
# demo_r_mwk_ts = Deaths by week and sex
z <- get_eurostat("demo_r_mwk_ts",  stringsAsFactors = FALSE) %>%
  mutate (time = as.character(time)) %>%
  mutate (year = as.numeric(substr(time, 1, 4)), 
          week = as.numeric(substr(time, 6,7)),
          sex = as.factor(sex),
          geo = as.factor(geo)) %>%
  select (sex, geo, year, week, value=values) %>%
  filter (sex == 'T')

Końcowy filter usuwa liczbę zgonów dla kobiet/mężczyzn osobno, bo chcemy tylko analizować ogółem czyli Total.

W pliku ecdc_countries_names_codes3166.csv są nazwy krajów (bo w bazie Eurostatu są tylko ISO-kody):

## country names
nn <- read.csv("ecdc_countries_names_codes3166.csv", sep = ';',  header=T, na.string="NA" )

Teraz liczę 5-letnie średnie dla okresu przed pandemią czyli dla lat 2015--2019 włącznie:

## mean weekly deaths 2015--2019
z0 <- z %>% filter ( year >= 2015 & year < 2020) %>% 
  group_by(geo, week) %>% 
  summarise (d0 = mean(value, na.rm=T))

Oraz tygodniowe liczby zgonów dla lat 2020--2021

## weekly deaths 2020--2021
z1 <- z %>% filter ( year > 2019 ) %>% 
  group_by(geo, year, week) %>% 
  summarise ( d1 = sum(value))

Łączę z0 z z1; liczę różnice między zgonami 2020--2021 a średnimi. Najpierw usuwam wiersze z NA w kolumnach zawierających zgony/średnie (drop_na). To w szczegolności usunie wiersze z tygodni (w roku 2021), dla których jeszcze nie ma danych.

## join z0 z1 and compute differences
zz <- left_join(z0, z1, by=c("week", "geo")) %>% 
  drop_na(d0,d1) %>%
  mutate (exp = (d1 - d0)/d0 * 100, exm = d1 - d0 )
zz <- left_join(zz, nn, by=c('geo'='id'))

Osobno liczę numer ostatniego raportowanego tygodnia w roku 2021. Pytanie o numer ostatniego w ogóle jest bez sensu, bo w 2020 będzie to ostatni tydzień. Ten fragment ciut mało robust jest bo w 2022 trzeba go będzie zmodyfikować

## if NA then the country stop reporting in 2020
zz.last.week <- zz %>% filter (year == 2021) %>% 
  group_by (geo) %>% 
  summarise (lw = last(week)) %>%
  drop_na(lw)
  
## najbardziej aktualny raportowany tydzień (różne kraje raportują różnie)
latestweek <- max (zz.last.week$lw)

Uwaga: niektóre kraje nie raportują w 2021, więc ich nie będzie w ramce zz.last.week z uwagi na filter.

Ramka zzt zawiera sumy różnic; bezwzględne (exm) oraz względne w procentach względem poziomu średniego 2015--2019 (exp):

### total exmort 
zzt <- zz %>% group_by(geo) %>%
   summarise(country=last(country), 
   exm = sum(exm),
   nm = sum(d0)) %>%
  mutate (exmp = exm/nm * 100)

Dodajemy informację z ramki zz.last.week; jeżeli w tej ramce nie ma kraju (bo nie raportował w 2021) to usuwamy go (drop_na(lw)). Wypierniczamy kraje raportujące więcej niż 6 tygodni po kraju (krajach), który przysłały raport najpóźniej. Np jak jakiś kraj-prymus dostarczył w 47 tygodniu to tylko kraje raportujące z tygodni 41 i później zostają a inne są usuwane (jako nieporównywalne)

zztt <- left_join(zzt, zz.last.week, by='geo') %>%
  drop_na(lw) %>%
  filter (lw >= latestweek - 6)
## ile zostało krajów (dla ciekawości):
countries.left <- nrow(zztt)

Wykres słupkowy względnej nadmiarowej umieralności (Total excess mortality by country) Trochę go komplikujemy bo słupek PL ma być wyróżniony (w tym celu tworzymy zmienną base a potem ręcznie ustawiamy legendę scale_fill_manual):

p6 <- zztt %>%  mutate( base=ifelse(geo=='PL', "1", "0")) %>% 
  ggplot(aes(x = reorder(country, exmp ), fill=as.factor(base))) +
  geom_bar(aes(y = exmp), stat="identity", alpha=.4 ) +
  geom_text(aes(label=sprintf("%.1f", exmp), y= exmp ), 
            vjust=0.25, hjust=1.25, size=2, color='black' ) +
  geom_text(aes(label=sprintf("%.0f", lw), y= 0 ), 
            vjust=0.25, hjust=-1.25, size=2, color='brown4' ) +
  xlab(label="") +
  ylab(label="") +
  ggtitle("Total Excessive deaths in Europe 2020--2021",
          subtitle ="Sum of (deaths_2020/2021 - average_2015--2019) / average_2015--2019 * 100%") +
  theme_nikw() +
  coord_flip() +
  scale_fill_manual( values = c( "1"="red", "0"="steelblue" ), guide = FALSE ) +
  labs(caption=x.note)

Wykres dynamiki. Ponieważ krajów jest dużo to usuwamy (pierwszy filter) `nieciekawe' (małe lub te które mają jakiś feler, np dane niepełne);

p7 <- zz %>% filter (! geo %in% c('AL', 'AM', 'IE', 'IS', 'UK', 'CY', 'GE', 'LI', 'ME')) %>%
  mutate(date= as.Date(sprintf("%i-%i-1", year, week), "%Y-%U-%u") ) %>%
  ##ggplot(aes(x = date, y =exp, group=geo, color=geo )) +
  ggplot(aes(x = date, y =exp)) +
  facet_wrap( ~country, scales = "fixed", ncol = 4, shrink = F) +
  geom_point(size=.4) +
  geom_smooth(method="loess", size=1, se = F, color="red",  span =.25) +
  geom_hline(yintercept = 50, color="firebrick", alpha=.2, size=1) +
  scale_y_continuous(breaks=seq(-100, 200, by=50)) +
  coord_cartesian(ylim = c(-100, 200)) +
  ggtitle("Excessive deaths in Europe 2020--2021",
          subtitle ="(deaths_2020/2021 - average_2015--2019) / average_2015--2019 * 100%") +
  theme_nikw() +
  xlab(label="%") +
  labs(caption=note)

Teraz wszystko na jednym ale ponieważ krajów jest za dużo (w sensie za mało byłoby kolorów żeby wykres był czytelny), to dzielimy je na dwie grupy: Polska i reszta. Każda grupa innym kolorem:

p8 <- zz %>% filter (! geo %in% c('AL', 'AM', 'IE', 'IS', 'UK', 'CY', 'GE', 'LI', 'ME')) %>%
  mutate( base=ifelse(geo=='PL', "1", "0")) %>% 
  mutate(date= as.Date(sprintf("%i-%i-1", year, week), "%Y-%U-%u") ) %>%
  ggplot(aes(x = date, y =exp,  color=as.factor(base ))) +
  ##ggplot(aes(x = date, y =exp)) +
  geom_point(size=.4) +
  geom_smooth(aes(x = date, y =exp, group=geo), method="loess", size=.3, se = F, span =.25) +
  geom_hline(yintercept = 50, color="firebrick", alpha=.2, size=1) +
  coord_cartesian(ylim = c(-100, 200)) +
  scale_color_manual( values = c( "1"="red", "0"="steelblue" ), guide = FALSE ) +
  theme_nikw() +
  ylab(label="%") +
  ggtitle("Excessive deaths in Europe 2020--2021 and Poland (red)",
          subtitle ="(deaths_2020/2021 - average_2015--2019) / average_2015--2019 * 100%") +
  scale_y_continuous(breaks=seq(-100, 200, by=50)) +
  labs(caption=note)

Wszystko działa (i koliduje jak mówił Bolec), bo działało na PC, ale był problem na raspberry, na którym to ostatecznie ma być uruchamiane z automatu co dwa tygodnie. Bo tam Debian starszy (Buster) i na arma. Zapodanie:

install.package("eurostat")

Skutkowało pobraniem wielu pakietów i błędem przy kompilacji czegoś co się nazywało s2. Kombinowałem na różne sposoby jak to skompilować bo stosownego pakietu w Debianie nie było. Ni-chu-chu.

W przypływie olśnienia obejrzałem zależności tego eurostatu, wśród których żadnego s2 nie było. Był za to sf, które z kolei rzekomo zależał od s2. Żeby było śmieszniej sf dało się zainstalować (mimo że s2 nie było -- stąd rzekomo):

apt install r-cran-sf/oldstable

Problem się rozwiązał...

Skrypt będzie się uruchamiał co dwa tygodnie. Generował rysunki i wysyłał je na twittera

Całość jest na githubie i to w dodatku jako plik Rmd, ale nie na moim koncie tylko na koncie Koła Naukowego SM w PSW w Kwidzynie, którego póki co jestem opiekunem i jedynym członkiem w jednym (por. github.com/knsm-psw/ES-mortality)

R w Debianie (wersja 10) jest w wersji 3.5:

> sessionInfo()
R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 10 (buster)

other attached packages:
[1] rmarkdown_1.17 ggplot2_3.2.1  ggthemes_4.2.

## https://cran.r-project.org/bin/linux/debian/#debian-buster-stable
## For a backport of R 4.0.1 to buster, please add
## deb http://cloud.r-project.org/bin/linux/debian buster-cran40/  
## /etc/apt/sources.list on your computer.
##
apt update

błąd...

## https://unix.stackexchange.com/questions/559977/not-able-to-install-latest-r-statistical-package-in-debian-stable
apt-key adv --keyserver keys.gnupg.net --recv-key FCAE2A0E115C3D8A
apt update
apt install -t buster-cran40 r-base

## Jako root
update.packages(lib.loc="/usr/local/lib/R/site-library", ask = FALSE,
  checkBuilt = TRUE, Ncpus = 16)

>  sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 10 (buster)

Wszystko działa do tej pory. Instaluję więcej pakietów:

install.packages("devtools")

Znowu błąd. Trzeba doinstalować pakiet libgit2-dev:

apt install libgit2-dev

install.packages("dplyr")
install.packages("tidyverse")
install.packages("ggthemes")
install.packages("bookdown")
install.packages("ggpubr")
install.packages("ggpubr")
install.packages("ggforce")
## or install.packages(c("ggpubr", "ggforce"))

Na sam koniec się okazuje, że pandoc (w wersji 2.2.1) stał się out-of-sync:

/usr/lib/rstudio/bin/pandoc/pandoc --version
pandoc 2.3.1
/usr/bin/pandoc --version
pandoc 2.2.1
## Po prostu
mv /usr/bin/pandoc /usr/bin/pandoc_2.2.1
ln -s /usr/lib/rstudio/bin/pandoc/pandoc /usr/bin/pandoc

Nie szło zainstalować (pakiet R) ggforce na raspberri, bo kompilacja kończyła się zwisem komputera. Coś mnie wreszcie tknęło, że problem może być związany z małą ilością pamięci

free
              total        used    
Mem:         948304      109108 
Swap:        102396       48680

swapon --show
NAME      TYPE SIZE  USED PRIO
/var/swap file 100M 47,6M   -2

Nie wiem kto wymyślił ten 100Mb swap, ale to bez sensu.

## tworzę drugi swap 
dd if=/dev/zero of=/var/swap.img bs=1024k count=1000
mkswap /tmp/swap.img
swapon /tmp/swap.img
free

Teraz ggforce się bezproblemowo instaluje.

## Usuwanie
swapoff -v /tmp/swap.img
sudo rm /tmp/swap.img

GUS się wychylił niespodziewanie z dużą paczką danych nt. zgonów. Dane są tygodniowe, w podziale na płcie, regiony w klasyfikacji NUTS oraz 5 letnie grupy wiekowe.

Dane są udostępnione w formacie XSLX w dość niepraktycznej z punktu widzenia przetwarzania strukturze (kolumny to tygodnie/wiersze to różne kategorie agregacji: płeć, wiek, region), który zamieniłem na CSV o następującej prostej 7 kolumnowej strukturze:

year;sex;week;date;age;geo;value

W miarę oczywiste jest, że year to rok, sex to płeć, week to numer tygodnia, date to data pierwszego dnia tygodnia (poniedziałek), geo to identyfikator obszaru a value liczba zgonów odpowiadająca kolumn 1--6. Ten plik jest podzielony na lata bo w całości zajmuje circa 200Mb. Umieściłem go tutaj.

Skrypt też w R wymodziłem co wizualizuje zgony wg grup wieku oraz województw. Ponieważ kombinacji płeć/wiek/region są setki, moje wykresy dotyczą zgonów ogółem/kobiet/mężczyzn w podziale na grupy wiekowe oraz ogółem w podziale na województwa. Każdy wykres zawiera dwa szeregi: liczbę zgonów w 2020 roku oraz średnią liczbę zgonów z lat 2015--2019. Ponadto jest wykres z jedną krzywą: procent liczony dla stosownych tygodni jako liczba zgonów w 2020 przez średnią 5 letnią z lat 2015--2019. Ten wykres występuje też w wariancie skróconym: tylko 6 ostatnich tygodni, co pozwala dodać do punktów wartości liczbowe (które nie zachodzą na siebie).

library("ggplot2")
library("dplyr")
library("scales")
library("ggthemes")
library("ggpubr")
library("tidyr")

picWd <- 12
spanV <- 0.5
GUS.url <- "https://stat.gov.pl/obszary-tematyczne/ludnosc/ludnosc/zgony-wedlug-tygodni,39,2.html"
NIKW.url <- "(c) NI-KW @ github.com/knsm-psw/GUS_mortality"
NIKW <- sprintf ("%s | %s", GUS, NIKW.url)

z <- read.csv("PL-mortality-2015.csv", sep = ';',  header=T, na.string="NA" )
lastO <- last(z$date)
lastT <- last(z$week)

nuts <- c('PL21', 'PL22', 'PL41', 'PL42', 'PL43', 'PL51', 'PL52', 'PL61', 'PL62',
'PL63', 'PL71', 'PL72', 'PL81', 'PL82', 'PL84', 'PL91', 'PL92')

### Ogółem
z00 <- z %>% filter ( sex == 'O'  & geo == 'PL' ) %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(age,week) %>% summarise( year = 't19',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
zz1 <- z1 %>% group_by(age,week) %>% summarise( year = 't20',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

farbe19 <- '#F8766D'
farbe20 <- '#00BFC4'

p1 <- ggplot(zz1, aes(x=week, y=vv, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~age, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"), values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg grup wiekowych (PL/Ogółem)", subtitle=sprintf("%s | ostatni tydzień: %s", NIKW, lastO) )

ggsave(plot=p1, "zgony_PL_by_age_O.png", width=picWd)

### M ###
z00 <- z %>% filter ( sex == 'M'  & geo == 'PL' ) %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(age,week) %>% summarise( year = 't19',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
zz1 <- z1 %>% group_by(age,week) %>% summarise( year = 't20',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

p2 <- ggplot(zz1, aes(x=week, y=vv, group=year, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~age, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 ##labs(caption=source) +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"),
   values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg grup wiekowych (PL/Mężczyźni)", subtitle=sprintf("%s | ostatni tydzień: %s", NIKW, lastO) )

ggsave(plot=p2, "zgony_PL_by_age_M.png", width=picWd)


### K #########################################
z00 <- z %>% filter ( sex == 'K'  & geo == 'PL' ) %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(age,week) %>% summarise( year = 't19',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
zz1 <- z1 %>% group_by(age,week) %>% summarise( year = 't20',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

p3 <- ggplot(zz1, aes(x=week, y=vv, group=year, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~age, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 ##labs(caption=source) +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"),
   values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg grup wiekowych (PL/Kobiety)", subtitle=sprintf("%s | ostatni tydzień: %s", NIKW, lastO) )

ggsave(plot=p3, "zgony_PL_by_age_K.png", width= picWd)

### ogółem wg województw #####################################
n <- read.csv("nuts.csv", sep = ';',  header=T, na.string="NA" )
## dodaj nazwy
z <- left_join(z, n, by='geo')

## wiek razem
z00 <- z %>% filter ( sex == 'O' & geo %in% nuts & age == 'OGÓŁEM') %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(name,week) %>%
 summarise( year = 't19', vv = mean(value, na.rm=TRUE)) %>% as.data.frame
 zz1 <- z1 %>% group_by(name,week) %>%
  summarise( year = 't20', vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

lastWeek <- last(zz1$week)
firstWeek <- lastWeek - 6

zz1 <- zz1 %>% filter ( week >= firstWeek  ) %>% as.data.frame
print(zz1)

p4 <- ggplot(zz1, aes(x=week, y=vv, group=year, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~name, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 ##labs(caption=source) +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"),
   values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg województw* (PL/Ogółem)", 
   subtitle=sprintf("*wg klasyfikacji NUTS stąd mazowieckie/stołeczne | %s | ostatni tydzień: %s", NIKW, lastO))

ggsave(plot=p4, "zgony_PL_by_woj_O.png", width=picWd)

## jako %% w średniej w poprzednich 5 lat

zz1 <- zz1 %>% spread(year, vv)

zz1$yy <- zz1$t20 / zz1$t19 * 100

p5 <- ggplot(zz1, aes(x=week, y=yy), color=farbe20) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4, color=farbe20) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~name, scales = "fixed") +
 xlab(label="nr tygodnia") +
 ylab(label="%") +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok 2020: ", labels = c("% 2020/(średnia 2015--2015)"),
   values = c("yy"=farbe20 )  ) +
 ggtitle("Zgony wg województw* (PL/Ogółem)", 
   subtitle=sprintf("*wg klasyfikacji NUTS stąd mazowieckie/stołeczne | %s | ostatni tydzień: %s", NIKW, lastO))

ggsave(plot=p5, "zgony_PL_by_woj_P.png", width=picWd)

zz1 <- zz1 %>% filter ( week >= firstWeek  ) %>% as.data.frame

p6 <- ggplot(zz1, aes(x=week, y=yy), color=farbe20) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4, color=farbe20) +
 geom_point(size=.4, alpha=.5) +
 geom_text(aes(label=sprintf("%.1f", yy)), vjust=-1.25, size=1.5) +
 facet_wrap( ~name, scales = "fixed") +
 xlab(label="nr tygodnia") +
 ylab(label="%") +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok 2020: ", labels = c("% 2020/(średnia 2015--2015)"),
   values = c("yy"=farbe20 )  ) +
   ggtitle(sprintf("Zgony wg województw* (PL/Ogółem) tygodnie: %i--%i (%i tydzień zaczyna się %s)",
     firstWeek, lastWeek, lastWeek, lastO), 
   subtitle=sprintf("*wg klasyfikacji NUTS stąd mazowieckie/stołeczne | %s", NIKW))

ggsave(plot=p6, "zgony_PL_by_woj_P6.png", width=picWd)

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".

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?

Reuters donosi, że While Russia has confirmed the world's fourth largest tally of coronavirus cases, it has a relatively low death toll from the associated disease, COVID-19[...] But data released by the Rosstat State Statistics Service on Sept. 4 show there were 57,800 excess deaths between May and July, the peak of the outbreak. The figure was calculated by comparing fatalities over those three months in 2020 with the average number of May-July deaths between 2015 and 2019. The excess total is more than three times greater than the official May-July COVID-19 death toll of 15,955. Co mnie zmotywowało do poszukania danych źródłowych.

Okazało się że dość prosto jest je pobrać ze strony Urzędu Statystycznego Federacji Rosyjskiej (https://eng.gks.ru/ a konkretnie https://showdata.gks.ru/finder/) Ponieważ okazało się to aż tak proste, to pobrałem nie tylko zgony ale także urodzenia. W okresie 2015--2020 (miesięcznie).

Dane są w formacie XSLX. Kolumny to miesiące wiersze poszczególne regiony i ogółem Federacja. Eksportuję do CSV używając LibreOffice i wybieram wiersz z danymi dla całej federacji. W rezultacie mam trzy wiersze: rok-miesiąc, urodzenia i zgony. Obracam moim super skryptem do transpozycji plików CSV:

month;urodzenia;zgony
2015-01-01;149269.99;174722.99
2015-02-01;144591.99;156737
2015-03-01;160974.99;175809.99
...

Trzy wykresy liniowe przedstawiające dynamikę zgonów, urodzin i przyrostu naturalnego (czyli różnicy między urodzinami a zgonami)

library("dplyr")
library("ggplot2")
library("scales")
spanV <- 0.5

## przyrost naturalny
d <- read.csv("UZRT.csv", sep = ';', header=T, na.string="NA");
d$diff <- d$urodzenia - d$zgony

pz <- ggplot(d, aes(x= as.Date(month), y=zgony)) + 
 geom_line(color="steelblue", size=.6, alpha=.5) +
 geom_point(color="steelblue", size=.8) +
 ## Trend 
 geom_smooth(method="loess", se=F, span=spanV, color="red3") +
 ## Skala na osi X-ów
 scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
 xlab(label="") +
 ylab(label="deaths") +
  geom_vline(xintercept = as.Date("2015-07-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2016-07-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2017-07-01"), alpha=.25, size=1)+
  geom_vline(xintercept = as.Date("2018-07-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2019-07-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2020-07-01"), alpha=.25, size=1) +
 labs(caption="https://showdata.gks.ru/finder/") +
 ggtitle("Russian Federation: Deaths (ths)",  
  subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: july", spanV)) 

pu <- ggplot(d, aes(x= as.Date(month), y=urodzenia)) + 
 geom_line(color="steelblue", size=.6, alpha=.5) +
 geom_point(color="steelblue", size=.8) +
 geom_smooth(method="loess", se=F, span=spanV, color="red3") +
 scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
 xlab(label="") +
 ylab(label="births") +
 labs(caption="https://showdata.gks.ru/finder/") +
  geom_vline(xintercept = as.Date("2015-08-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2016-08-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2017-08-01"), alpha=.25, size=1)+
  geom_vline(xintercept = as.Date("2018-08-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2019-08-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2020-08-01"), alpha=.25, size=1) +
 ggtitle("Russian Federation: Births (ths)",
  subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: august", spanV)) 

pp <- ggplot(d, aes(x= as.Date(month), y=diff)) + 
 geom_line(color="steelblue", size=.6, alpha=.5) +
 geom_point(color="steelblue", size=.8) +
 geom_smooth(method="loess", se=F, span=spanV, color="red3") +
 scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
  geom_vline(xintercept = as.Date("2015-05-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2016-05-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2017-05-01"), alpha=.25, size=1)+
  geom_vline(xintercept = as.Date("2018-05-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2019-05-01"), alpha=.25, size=1) +
  geom_vline(xintercept = as.Date("2020-05-01"), alpha=.25, size=1) +
 xlab(label="") +
 ylab(label="balance") +
 labs(caption="https://showdata.gks.ru/finder/") +
 ggtitle("Russian Federation: natural balance (ths)", 
  subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: may", spanV))

Wykres liniowy dla miesięcy jednoimiennych (MJ):

library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
options(dplyr.print_max = 1e9)
size0 <- 1.2
size1 <- 0.6
size2 <- 0.6

tt <- read.csv("UZRT.csv", sep = ';',  header=T, na.string="NA");
tt$rok <- substr(tt$month,  1, 4)

tt$yyyymmdd <- as.Date(sprintf ("%s-%s", substr(tt$month,  6, 7),
substr(tt$month,  9, 10) ), format="%m-%d")

## Obliczenie przyrostu naturalnego
tt$diff <- tt$urodzenia - tt$zgony

## Obliczenie 
tt.yy2020 <- tt %>% filter(rok==2020) %>% as.data.frame
tt.yy2019 <- tt %>% filter(rok==2019) %>% as.data.frame
tt.yy2018 <- tt %>% filter(rok==2018) %>% as.data.frame
tt.yy2017 <- tt %>% filter(rok==2017) %>% as.data.frame
tt.yy2016 <- tt %>% filter(rok==2016) %>% as.data.frame
tt.yy2015 <- tt %>% filter(rok==2015) %>% as.data.frame

pf <- ggplot() +
    ggtitle("Russian Federation: births 2015--2020") +
    geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2020'), alpha=.5, size=size0) +
    geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2020'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2019'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2019'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2018'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2018'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2017'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2017'), alpha=.5, size=size1) +
    ###
    geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2016'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2016'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2015'), alpha=.5, size=size2) +
    geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = urodzenia,  colour = '2015'), alpha=.5, size=size1) +
    scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
    labs(caption="https://showdata.gks.ru/finder/") +
    ylab(label="ths") +
    xlab(label="") +
    labs(colour = "") +
    theme(legend.position="top") +
    theme(legend.text=element_text(size=10));

qf <- ggplot() +
    ggtitle("Russian Federation: deaths 2015-2020") +
    geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2020'), alpha=.5, size=size0) +
    geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2020'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2019'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2019'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2018'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2018'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2017'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2017'), alpha=.5, size=size1) +
    ###
    geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2016'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2016'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2015'), alpha=.5, size=size2) +
    geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = zgony,  colour = '2015'), alpha=.5, size=size2) +
    scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
    labs(caption="https://showdata.gks.ru/finder/") +
    ylab(label="ths") +
    xlab(label="") +
    labs(colour = "") +
    theme(legend.position="top") +
    theme(legend.text=element_text(size=10));

qf <- ggplot() +
    ggtitle("Russian Federation: natural balance 2015-2020") +
    geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = diff,  colour = '2020'), alpha=.5, size=size0) +
    geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = diff,  colour = '2020'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = diff,  colour = '2019'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = diff,  colour = '2019'), alpha=.5, size=size0) +
    geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = diff,  colour = '2018'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = diff,  colour = '2018'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = diff,  colour = '2017'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = diff,  colour = '2017'), alpha=.5, size=size1) +
    ###
    geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = diff,  colour = '2016'), alpha=.5, size=size1) +
    geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = diff,  colour = '2016'), alpha=.5, size=size1) +
    geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = diff,  colour = '2015'), alpha=.5, size=size2) +
    geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = diff,  colour = '2015'), alpha=.5, size=size2) +
    scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
    labs(caption="https://showdata.gks.ru/finder/") +
    ylab(label="ths") +
    xlab(label="") +
    labs(colour = "") +
    theme(legend.position="top") +
    theme(legend.text=element_text(size=10));

Wreszcie wykres kombinowany liniowo-słupkowy dla miesięcy jednoimiennych (MJ). Linie to 1) średnia liczby zgonów w latach 2015--2020, 2) średnia liczba zgonów w latach 2015--2020 plus/minus odchylenie standardowe liczby zgonów (dla MJ); 3) liczba zgonów w roku 2020. Wykres słupkowy to różnica między liczbą zgonów w roku 2020 a maksymalną liczbą zgonów w latach 2015-2019 ORAZ liczba zgonów z powodu COVID19 (z bazy ECDC). Ten skrypt korzysta z danych z pliku covid_ru.csv, który powstał przez zagregowanie danych dziennych dostępnych ze strony ECDC (dla Federacji Rosyjskiej) oraz policzenie średnich/odchyleń dla jednoimiennych miesięcy z danych w pliku UZRT.csv

library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
options(dplyr.print_max = 1e9)

#
tt <- read.csv("covid_ru.csv", sep = ';',  header=T, na.string="NA");
tt$diff <- tt$urodzenia - tt$zgony

cols <- c("mean19"="navyblue","mean+sd19"="#3591d1", "mean-sd19"="#3591d1",
  "deaths20"='brown4', "diff19"="red", "c19d"='cyan')

pf <- ggplot(data=tt) +
    ggtitle("Russian Federation: deaths 2015--2020", 
             subtitle='mean19/sd19: mean/sd 2015--19; diff19: 2020 - max(2015-019); deaths20: deaths 2020; c19d: C19 deaths') +
    geom_line(mapping = aes(x=yyyymmdd, y = mean,  colour = 'mean19'), alpha=.5, size=.8) +
    geom_line(mapping = aes(x=yyyymmdd, y = mean + sd,  colour = 'mean+sd19'), alpha=.5, size=.8) +
    geom_line(mapping = aes(x=yyyymmdd, y = mean - sd,  colour = 'mean-sd19'), alpha=.5, size=.8) +
    geom_text(mapping = aes(x=yyyymmdd, label=difflabel, y=diff), vjust=-1.0, size=3 ) +
    ##geom_line(mapping = aes(x=yyyymmdd, y = c19d,  colour = 'c19d'), alpha=.5, size=1.2) +
    labs(caption = "Source: https://showdata.gks.ru/finder/; https://www.ecdc.europa.eu/en/covid-19-pandemic") +
    ##
    geom_bar(mapping = aes(x=yyyymmdd, y = diff, fill = 'diff19' ), color='brown', stat="identity", alpha=.25) +
    geom_bar(mapping = aes(x=yyyymmdd, y = c19d, fill = 'c19d' ), color ='navyblue', stat="identity", alpha=.25) +
    ###
    geom_line(mapping = aes(x=yyyymmdd, y = deaths,  colour = 'deaths20'), alpha=.5, size=1.2) +
    scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
    ylab(label="ths") +
    xlab(label="") +
    ##labs(colour = "Legend") +
    scale_colour_manual(name="Lines: ", values=cols) +
    scale_fill_manual(name="Bars: ", values=cols) +
    theme(legend.position="right") +
    theme(legend.text=element_text(size=12));

Skrypt pomocniczy liczący średnie 2015--2019 dla miesięcy jednoimiennych itp rzeczy

library(dplyr)
library(tidyr)
## Agreguje zgodny wg miesięcy jednoimiennych (MJ) dla lat 2015-19
## drukuje średnią/sd/max dla MJ
## Drukuje z2020 - max(z2015--2019)

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

tt$rok <- substr(tt$month,  1, 4)
tt$yyyymmdd <- as.Date(sprintf ("%s-%s", substr(tt$month,  6, 7), substr(tt$month,  9, 10) ), format="%m-%d")

tt$diff <- tt$urodzenia - tt$zgony

tt.yy2020 <- tt %>% filter(rok==2020) %>% as.data.frame
tt.yy2019 <- tt %>% filter(rok==2019) %>% as.data.frame
tt.yy2018 <- tt %>% filter(rok==2018) %>% as.data.frame
tt.yy2017 <- tt %>% filter(rok==2017) %>% as.data.frame
tt.yy2016 <- tt %>% filter(rok==2016) %>% as.data.frame
tt.yy2015 <- tt %>% filter(rok==2015) %>% as.data.frame

## max zgony
tt.yy2020$max.so.far <-  pmax (tt.yy2019$zgony, tt.yy2018$zgony,
  tt.yy2017$zgony, tt.yy2016$zgony, tt.yy2015$zgony)

## średnia
tt.yy2020$mean <-  (tt.yy2019$zgony + tt.yy2018$zgony
  + tt.yy2017$zgony + tt.yy2016$zgony + tt.yy2015$zgony)/5

## odchylenie std:
tt.yy2020$sqmean <-  (tt.yy2019$zgony^2 + tt.yy2018$zgony^2 +
   tt.yy2017$zgony^2 + tt.yy2016$zgony^2 + tt.yy2015$zgony^2)/5
tt.yy2020$sq <- sqrt(tt.yy2020$sqmean - tt.yy2020$mean^2)

tt.yy2020$diff20 <- tt.yy2020$zgony - tt.yy2020$max.so.far

sprintf ("%s;%.2f;%.2f;%.2f;%.2f", 
    tt.yy2020$yyyymmdd, tt.yy2020$diff20, tt.yy2020$sqmean, tt.yy2020$sq, tt.yy2020$max.so.far )

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:

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

Trzeba coś robić w czasie kwarantanny

## https://b-rodrigues.github.io/modern_R/
## https://gist.github.com/imartinezl/2dc230f33604d5fb729fa139535cd0b3
library("eurostat")
library("dplyr")
library("ggplot2")
library("ggpubr")
## 
options(scipen=1000000)
dformat <- "%Y-%m-%d"

eu28 <- c("AT", "BE", "BG", "HR", "CY", "CZ", "DK",
   "EE", "FI", "FR", "DE", "EL", "HU", "IE", 
   "IT", "LT", "LU", "LV", "MT", "NL", "PL", 
   "PT", "RO", "SK", "SI", "ES", "SE")
eu6 <- c("DE", "FR", "IT", "ES", "PL")

### Demo_mor/ Mortality monthly ### ### ###
dm <- get_eurostat(id="demo_mmonth", time_format = "num");
dm$date <- sprintf ("%s-%s-01", dm$time, substr(dm$month, 2, 3))
str(dm)

## There are 12 moths + TOTAL + UNKN
dm_month <- levels(dm$month)
dm_month

## Only new data
dm28  <- dm %>% filter (geo %in% eu28 & as.Date(date) > "1999-12-31")
str(dm28)
levels(dm28$geo) 

## Limit to DE/FR/IT/ES/PL:
dm6  <- dm28 %>% filter (geo %in% eu6)
str(dm6)
levels(dm6$geo) 

pd1 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
 geom_line(aes(group = geo, color = geo), size=.4) +
 xlab(label="") +
 ##scale_x_date(date_breaks = "3 months", date_labels = "%y%m") +
 scale_x_date(date_breaks = "6 months",
   date_labels = "%m\n%y", position="bottom") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Deaths", subtitle="https://ec.europa.eu/eurostat/data/database (demo_mmonth)")

## Newer data
dm6  <- dm6 %>% filter (as.Date(date) > "2009-12-31")

pd2 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
 geom_line(aes(group = geo, color = geo), size=.4) +
 xlab(label="") +
 scale_x_date(date_breaks = "3 months", date_labels = "%m\n%y", position="bottom") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Deaths", subtitle="https://ec.europa.eu/eurostat/data/database (demo_mmonth)")

ggsave(plot=pd1, file="mort_eu_L.png", width=12)
ggsave(plot=pd2, file="mort_eu_S.png", width=12)

## Live births (demo_fmonth) ### ### ###

dm <- get_eurostat(id="demo_fmonth", time_format = "num");
dm$date <- sprintf ("%s-%s-01", dm$time, substr(dm$month, 2, 3))
str(dm)

## There are 12 moths + TOTAL + UNKN
dm_month <- levels(dm$month)
dm_month

dm28  <- dm %>% filter (geo %in% eu28 & as.Date(date) > "1999-12-31")
str(dm28)
levels(dm28$geo) 

dm6  <- dm28 %>% filter (geo %in% eu6)
str(dm6)
levels(dm6$geo) 

pd1 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
 geom_line(aes(group = geo, color = geo), size=.4) +
 xlab(label="") +
 ##scale_x_date(date_breaks = "3 months", date_labels = "%y%m") +
 scale_x_date(date_breaks = "6 months", date_labels = "%m\n%y", position="bottom") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Births", subtitle="https://ec.europa.eu/eurostat/data/database (demo_fmonth)")

##
dm6  <- dm6 %>% filter (as.Date(date) > "2009-12-31")

pd2 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
 geom_line(aes(group = geo, color = geo), size=.4) +
 xlab(label="") +
 scale_x_date(date_breaks = "3 months", date_labels = "%m\n%y", position="bottom") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Births", subtitle="https://ec.europa.eu/eurostat/data/database (demo_fmonth)")

ggsave(plot=pd1, file="birt_eu_L.png", width=12)
ggsave(plot=pd2, file="birt_eu_S.png", width=12)

## Population (only) yearly ### ### ##
## Population change - Demographic balance and crude rates at national level (demo_gind)
dp <- get_eurostat(id="demo_gind", time_format = "num");
dp$date <- sprintf ("%s-01-01", dp$time)
str(dp)
dp_indic_dic <-  get_eurostat_dic("indic_de")

dp_indic_dic
dp28  <- dp %>% filter (geo %in% eu28 & time > 1999 & indic_de == "JAN")

str(dp28)
dp6  <- dp28 %>% filter (geo %in% eu6)

pdp1 <- ggplot(dp6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
        geom_line(aes(group = geo, color = geo), size=.4) +
        xlab(label="") +
        ##scale_x_date(date_breaks = "3 months", date_labels = "%y%m") +
        ##scale_x_date(date_breaks = "6 months", date_labels = "%m\n%y", position="bottom") +
        theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
        ggtitle("Population", subtitle="https://ec.europa.eu/eurostat/data/database (demo_fmonth)")

ggsave(plot=pdp1, file="pdp1", width=12)

Na stronie http://policja.pl/pol/form/1,Informacja-dzienna.html udostępniane są dzienne dane dotyczące liczby interwencji, zatrzymanych na gorącym uczynku, zatrzymanych poszukiwanych, pijanych kierujących, wypadków, zabitych w wypadkach, rannych w wypadkach.

Ściągam wszystkie dane:

#!/bin/bash

rm pp.html

for ((i=0; i <= 274; i++)) do 
  if [ ! -f ${i}.html ] ; then
    curl -o ${i}.html "http://policja.pl/pol/form/1,Informacja-dzienna.html?page=${i}" ; 
    grep 'data-label' ${i}.html >> pp.html
    sleep 6
  else 
    grep 'data-label' ${i}.html >> pp.html
    echo done
  fi

done

zamieniam prostymi skryptami na plik CSV, który ma następującą strukturę:

data;interwencje;zng;zp;znk;wypadki;zabici;ranni
2008-12-01;NA;873;344;447;135;1;1

okazuje się że liczba interwencji jest podawana od roku 2018, wcześniej nie była. Nic to wstawiamy NA.

Na przyszłość dane będą aktualizowane w ten sposób, że codziennie (przez odpowiedni wpis w pliku crontab) będzie pobierany plik http://policja.pl/pol/form/1,Informacja-dzienna.html:

#!/usr/bin/perl
use LWP::Simple;

$PP="http://policja.pl/pol/form/1,Informacja-dzienna.html";
$PPBase="pp.csv";

$content = get("$PP");

$content =~ s/\r//g; # dla pewności usuń

@content = split (/\n/, $content);

foreach (@content) { chomp();
  unless ($_ =~ m/data-label=/ ) { next }

  if ($_ =~ m/Data statystyki/ ) { $d = clean($_); }
  elsif ($_ =~ m/Interwencje/ )  { $i = clean($_); }
  elsif ($_ =~ m/Zatrzymani na g/ ) { $zg = clean($_); }
  elsif ($_ =~ m/Zatrzymani p/ ) { $zp = clean($_); }
  elsif ($_ =~ m/Zatrzymani n/ ) { $zn = clean($_); }
  elsif ($_ =~ m/Wypadki d/ ) { $w = clean($_);  }
  elsif ($_ =~ m/Zabici/ )  { $z = clean($_);  }
  elsif ($_ =~ m/Ranni/ ) { $r = clean($_);
    $l = "$d;$i;$zg;$zp;$zn;$w;$z;$r";
    $last_line = "$l"; $last_date = "$d";
    ## pierwszy wpis powinien zawierać dane dotyczące kolejnego dnia
    ## więc po pobraniu pierwszego można zakończyć
    last;
 }
}

### read the database
open (PP, "<$PPBase") || die "cannot open $PPBase/r!\n" ;

while (<PP>) { chomp(); $line = $_; @tmp = split /;/, $line; }

close(PP);

### append the database (if new record)
open (PP, ">>$PPBase") || die "cannot open $PPBase/w!\n" ;

unless ("$tmp[0]" eq "$last_date") { print PP "$last_line\n" }
else {print STDERR "nic nowego nie widzę!\n"}

close(PP);

sub clean  {
 my $s = shift;
 $s =~ s/<[^<>]*>//g;
 $s =~ s/[ \t]//g;

 return ($s);
}

Zaktualizowana baza jest wysyłana na githuba. Tutaj jest: https://github.com/hrpunio/Nafisa/tree/master/PP

Agregacja do danych tygodniowych okazała się nietrywialna

Niektóra lata zaczynają się od tygodnia numer 0 a inne od 1. Okazuje się, że tak ma być (https://en.wikipedia.org/wiki/ISO_week_date#First_week):

If 1 January is on a Monday, Tuesday, Wednesday or Thursday, it is in W01. If it is on a Friday, it is part of W53 of the previous year. If it is on a Saturday, it is part of the last week of the previous year which is numbered W52 in a common year and W53 in a leap year. If it is on a Sunday, it is part of W52 of the previous year.

Nie bawię się w subtelności tylko tygodnie o numerze zero dodaję do tygodnia z poprzedniego roku.

Sprawdzam czy jest OK i się okazuje że niektóre tygodnie mają 8 dni. W plikach html są błędy:

Błędne daty 2019-10-30 winno być 2019-09-30; podobnie błędne 2019-03-28 (winno być 2019-02-28), 2018-11-01 (2018-12-01), 2018-12-01 (2017-12-01), 2016-04-30 (2016-03-30), 2009-08-31 (2009-07-31). Powtórzone daty: 2016-03-10, 2010-07-25, 2010-01-10 (zdublowane/różne/arbitralnie usuwamy drugi) Ponadto brak danych z następujących dni: 2015-12-04--2015-12-07, 2015-04-17--2015-04-20, 2014-10-02--2014-10-05, 2014-01-23 i jeszcze paru innych (nie chcialo mi się poprawiać starych.)

Teraz jest OK, plik ppw.csv ma nast strukturę:

rok;nrt;interwencje;in;zng;zngn;zp;zpn;znk;znkn;wypadki;wn;zabici;zn;ranni;rn;d1;d7 coś co się kończy na `n' to liczba tego co jest w kolumnie poprzedniej, np zn to liczba dni tygodnia dla kolumny zabici. Generalnie kolumny kończące się na `n' zawierają 7 :-) Kolumna d1 to pierwszy dzień tygodnia a kolumna d7 ostatni.

maxY <- max (d$zabici)
pz <- ggplot(d, aes(x= as.factor(nrt), y=zabici )) + 
 geom_bar(fill="steelblue", stat="identity")  +
 xlab(label="") +
 ggtitle("Wypadki/zabici (Polska/2020)", subtitle="policja.pl/pol/form/1,Informacja-dzienna.html") 

W sumie agregacja jest niepotrzebna, bo można ją zrobić na poziomie R używając funkcji stat_summary:

pw <- ggplot(d, aes(x= week, y=wypadki)) + 
 stat_summary(fun.y = sum, geom="bar", fill="steelblue") +
 scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
 xlab(label="") +
 #theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Wypadki (Polska/2018--2020)", subtitle="policja.pl/pol/form/1,Informacja-dzienna.html") 

albo najpierw agregując dane a potem wykreślając wykres szeregu zagregowanego. Drugi sposób pozwala na przykład na dodanie linii oznaczających poziomy zagregowanego zjawiska/etykiety słupków w sposób `inteligentny'. Dodajemy etykiety (z numerem tygodnia) tylko dla słupków poniżej/powyżej Q1/Q3:

## agregowanie do danych tygodniowych kolumn ranni, zabici, wypadki
dw <- d %>% group_by ( YrWeek) %>% summarise_at ( vars(ranni,zabici,wypadki), sum )

## Obliczanie mediany i kwartyli
median.zw <- median(dw$zabici)
q1.zw <- quantile(dw$zabici, 1/4) 
q3.zw <- quantile(dw$zabici, 3/4) 

## YrWeekZZ to numer tygodnia, dla tygodni w których liczba wypadków jest typowa
## numer jest pusty (żeby nie był drukowany); taki trick może jest lepszy
dw$YrWeekZZ <- substr(dw$YrWeek,4,5)
dw$YrWeekZZ[ (dw$zabici > q1.zw) & (dw$zabici < q3.zw) ] <- ""

pz.2 <- ggplot(dw, aes(x= YrWeek, y=zabici)) + 
 geom_bar(stat="identity", fill = "steelblue") +
 geom_text(data=dw, aes(label=sprintf("%s", YrWeekZZ), x=YrWeek, y= zabici), vjust=-0.9, size=3 ) +
 geom_hline(yintercept=median.zw, linetype="solid", color = "violet", size=.4) +
 geom_hline(yintercept=q3.zw, linetype="solid", color = "red", size=.4) +
 geom_hline(yintercept=q1.zw, linetype="solid", color = "red", size=.4) +
 xlab(label="rok/tydzień") +
 ylab(label="zabici") +
 scale_x_discrete(breaks=c("18/01", "18/10", "18/20",  "18/30", "18/40",
          "19/01", "19/10", "19/20",  "19/30", "19/40", "20/01", "20/10"))  +
          #  labels = c("/18/01", "18/10", "18/20", "")) ## tutaj niepotrzebne
 ggtitle("Wypadki/zabici (Polska/2018--2020)", 
  subtitle="Linie poziomie: q1/me/q3 (źródło: policja.pl/pol/form/1,Informacja-dzienna.html)") 

29 czerwca był upał i pierwszy raz w życiu zobaczyłem stację MEVO literalnie zawaloną rowerami. Była 13:00. Po dojechaniu do domu sprawdziłem, że rowerów było tam aż 24. To mnie zainspirowało do sprawdzenia jak wygląda koncentracja rowerów na stacjach. Na szybko zmajstrowałem skrypt wypisujący ile jest rowerów na stacjach (rowery), udział w całości zaparkowanych w danym momencie (udzial) oraz udział w całości zaparkowanych w danym momencie w tym konkretnym mieście (ludzial): Dla 29 czerwca 2019, godz 13:00 okazało się że:

stacja;wspolrzedne;miasto;rowery;udzial;ludzial
S11358;18.57196808,54.40258423;Gdańsk;55;6.782;16.566
S11069;18.59038998,54.42794681;Gdańsk;24;2.959;7.229
S10100;18.57092997,54.44534256;Sopot;24;2.959;35.294
S11357;18.63640600,54.38637008;Gdańsk;18;2.219;5.422
S12007;18.53351997,54.49663118;Gdynia;16;1.973;7.843
S11186;18.57639103,54.39859534;Gdańsk;15;1.850;4.518
S10121;18.56255831,54.45392821;Sopot;13;1.603;19.118
S12126;18.47267507,54.54728035;Gdynia;13;1.603;6.373
S12053;18.54832346,54.51633152;Gdynia;13;1.603;6.373
S12054;18.54687652,54.51960321;Gdynia;12;1.480;5.882
S12033;18.56357898,54.48005340;Gdynia;10;1.233;4.902

Czyli na stacji 11358 było 55 rowerów co stanowiło 6,782% wszystkich zaparkowanych w systemie MEVO albo 16,566% zaparkowanych w Gdańsku. Dla Sopotu było nawet jeszcze lepiej bo na stacjach 10100/10121 było 24+13 (37 rowerów) ale było to 35,294 + 19.118, tj. prawie 55% wszystkich zaparkowanych w Sopocie (o godzinie 13:00). Dokładnie było 68 rowerów w 28 miejscach wtedy (27 stacji i jeden luźny bajk). Na 11 stacjach nie było nic a na 10 jeden rower.

Jest taka prosta miara koncentracji, co się nazywa w języku HH-Index, albo po polsku Wskaźnik Herfindahla-Hirschmana. Jest on nawet stosowany w USA do mierzenia koncentracji na rynku. Formuła jest banalnie prosta: dla $N$ wartości $x_i$ ($i=1...N$) sumujących się do 100 (czyli udziałów w całości), HHI liczy się jako: $\sum_{i=1}^N x_i^2$. Łatwo sprawdzić że HHI < 10000. Interpretacja jest taka, że HHI < 1000 wskazuje na słabą koncentrację, 1000 < HHI < 1800 umiarkowaną, a wartości większe od 1800 na dużą. BTW HHI da Sopotu o 13:00 (29/7/2019) wynosiło około 1850...

No to ja policzyłem HHI dla MEVO. Udziały były chwilowe, tj. $x_i = r_i/r_t$, dzie $r_t$ łączna liczba zaparkowanych rowerów na wszystkich stacjach w mieście $M$ (w danym momencie); no a $x_i$ to oczywiście liczba rowerów na stacji $i$. Potem uśredniłem, tj. wszystkie $HHI_i$ z godziny $h$ zsumowałem i podzieliłem przez liczbę pomiarów w tej godzinie (zwykle przez 30, bo pomiar jest co 2 minuty). Policzyłem oddzielnie dla GD/GA/Sopot/Tczew dla dni pracujących oraz dla świąt, sobót i niedziel osobno...

Wyniki dla maj--lipiec na wykresie (obok). Ciekawostkowo koncentracja w GD jest zaskakująco inna niż w GA. Teoretycznie im więcej stacji tym wartość HHI powinna być mniejsza, a tak nie jest: w GD wartość HHI wynosi w szczycie około 500, a w GA tylko 300. Szczyt wypada tak 9--11 zresztą. Można sobie wyobrazić, że użytkownicy jadą do pracy i zostawiając rowery przez biurami ogołacają stacje poza centrum a zapełniają te w rodzaju stacji o numerze 11358 (Gdańsk/Oliva Business Centre). W niedziele i święta do pracy nie jeżdżą to koncentracja jest mniejsza. Ma sens, ale nie w GA gdzie akurat w święta jest większa, wprawdzie chwilowo (w znaczeniu, że szybko rośnie ale potem równie szybko spada), ale jednak. Jakby w GA masowo jeździli gdzieś koło południa, a potem wracali z powrotem prawie że od razu... Mniejsza liczba stacji/rowerów w GA niż GD może powodować że wartość HHI ,,łatwiej'' rośnie...

W Sopocie i Tczewie jest znacząco mniej rowerów niż w GD/GA więc nic dziwnego że wartość HHI jest też dużo większa. BTW przeciętne dzienne wartości HHI są następujące (GD/GA/Sopot/Tczew): 170/222/1255/1143 (pon-piątek) oraz 93/225/1169/1200 (niedziele-soboty-święta). W Sopocie (poniedziałek--piątek) zmiany HHI są jeszcze inne niż w GD/GA. Amplituda jest mniejsza, a jedyny wyraźny dół jest rano około 5 a nie w godzinach 21--5 jak na przykład w GD. W niedziele i święta jest tradycyjnie jeden szczyt koło południa, a oprócz tego około 19--20 oraz 2--3 rano.

Skrypty i plik CSV z danymi jest tradycyjnie w archiwum GitHub

Na koniec przypomnienie, że prezes obiecał na 18.08 +4000 rowerów w systemie. Na koniec lipca wykazywane jest póki co: +1500 (niektóre w Warszawie albo w Cedrach Wlk.). Z tej puli bajków (numerów bajków?) codziennie pojawia się w pliku locations.js +1350, jeździ zaś mniej, około 1200 (reszta stoi). Na 100 procent wszyscy już wiedzą, że ni-chu-chu nie będzie 4000 rowerów nawet na 31 września, ale media zachowują w tej sprawie zgodne milczenie. Nikt nikogo nie pyta, a zwłaszcza prezesa. Nie ma sprawy...

Powyborczo. Primo: szkoda, że PiSom te kamery nie wyszły byłoby się z czego pośmiać oglądając debili ze świeczkami (innego pożytku z zainstalowania -- przy założeniu #1kamera na jedną komisję -- 27tys kamer nie widzę).

Ale do rzeczy: dane pobrane z PKW (na Wikipedii za 2014 mają dokładnie takie same, za 2018 nie sprawdzałem)

require(ggplot2)

#d <- read.csv("mandaty.csv", sep = ';',  header=T, na.string="NA");
# Albo po prostu bo danych mało
# https://www.datamentor.io/r-programming/data-frame/
x <- data.frame("komitet" = c("PIS", "PO", "PSL", "SLD", "INNI"),
   "y2018" = c(254,194,70,15,19),
   "y2014" = c(171,179,157,28,20) )

# różnica w liczbie uzyskanych mandatów
d$diff <- d$y2018 - d$y2014 

ggplot(d, aes(x= komitet, y=diff, fill=komitet)) +
  geom_bar(stat="identity") +
  scale_fill_manual("legend",
    values = c("PIS" = "#421C52", "PO" = "blue",
    "PSL" = "green", "SLD" = "red", "INNI" = "pink")) +
    geom_text(aes(x=komitet, y=diff, label=diff),
    hjust=0, vjust=-0.25, size=3.5) +
ggtitle ("Mandaty sejmików wojewódzkich 2018--2014 (zmiana)")

BTW, jeżeli protokoły komisje obwodowe wysłały (zapewne elektronicznie) do PKW góra w poniedziałek (w mojej już poniedziałek-rano okleili kopią drzwi), to co niezawisłe Hermelińskie robiły we wtorek i środę? Niestety tego prostego pytania żaden z tzw. dziennikarzy (aka specjalistów od pierdołowatych njusów czyli #pierdokontentu) nie zadał.

A mnie ono ciekawi.

BTW2: ten wpis jest 500 w blogu. Wychodzi jakieś 45/rok średnio (z tendencją spadkową).

Analiza eksploracyjna wyborów wójtów/burmistrzów/prezydentów. W PL wybiera się radnych w wyborach do rad powiatów/rad gmin (oba ciała IMO zbędne), radnych sejmików wojewódzkich oraz uwaga: wójtów/burmistrzów/prezydentów (WBP) na poziomie gmin. O ile wybory sejmików kierują się tym samym mechanizmem co wybory sejmowe to wybory WBP są większościowe -- każdy może wystartować i wygrać. Do tego taki WBP ma dużą władzę więc warto być WBP. Takich wyborów w PL jest 2477 -- tyle ile gmin. W zależności od statusu gminy w jednych wybiera się wójta a w innych prezydenta czy burmistrza. Mówiąc konkretnie wójtów jest 1547, burmistrzów 823 a prezydentów 107. Poniższa tabela zestawia dane dotyczące kandydatów w wyborach 2018/2014/2010

  Rok      N   1KN   1K%   2KN   2K%   >4N    >4%     śr
  ------------------------------------------------------
  2018  6965   329  13,30  849  34,30  262  10,58   2,81
  2014  8019   245   9,90  666  26,90  471  19,00   3,25
  2010  7776   303  12,20  683  27.57  430  17,36   3,14

Jak na moje to kandydatów za dużo to się nie zgłasza, do tego (w tym roku/w tych wyborach) w 13,30% gmin jest jeden, a w 34,30% dwóch (co daje co najwyżej dwóch w prawie połowie wyborów WBP). Do tego tendencja jest jakby nie w tę stronę co trzeba: mniej kandydatów ogółem, więcej gmin z małą liczbą kandydatów, mniej gmin z dużą liczbą kandydatów. Można podsumować że demokracja na lokalnym poziomie słabnie...

Ilustruje to wykres krzywych gęstości liczby kandydatów na urząd WBP (dla każdego roku oddzielna krzywa).

## ramka g ma następującą strukturę: razem;teryt;rok
g$r <- as.factor(g$rok)
p <- ggplot(g, aes(x=razem, color=r)) + geom_density() +
labs(title="Krzywa gęstości liczby kandydatów na urząd wójta/burmistrza/prezydenta",
x="Liczba kandydatów", 
y = "Gęstość", color="Rok")

Dane są tutaj

Dobrnąłem w końcu do finału pobierając ostatecznie ze strony PKW dane dotyczące siedmiu wyborów, które odbyły się w latach: 2015, 2014 (samorządowe), 2011, 2010 (samorządowe), 2007, 2006 (samorządowe), 2005.

Wyniki wcześniejszych wyborów nie są już dostępne na poziomie komisji obwodowych (a przynajmniej ja nie potrafię takowych odszukać). Protokoły z wyborów z 2006 roku też nie były dostępne, ale udało się je w części odtworzyć ze stron z wynikami kandydatów (zawierającymi liczbę głosów oddanych na kandydata w poszczególnych komisjach obwodowych).

Dla każdych wyborów wykreśliłem histogram poparcia dla mainstreamowych partii: PSL, PO, PiS oraz SLD. Zgodnie z oczekiwaniami rozkłady poparcia są jednomodalne, prawostronnie symetryczne, ale z dwoma wyjątkami: rozkład poparcia dla PO jest bimodalny i ta tendencja wydaje się stała. Rozkład poparcia dla PSL z roku 2014 (cud nad urną) różni się -- na zasadzie znajdź element niepasujący do pozostałych -- od sześciu pozostałych rozkładów poparcia dla tej partii (czemu to już inna historia).

Dane są tutaj

Rozkład odsetka głosów nieważnych (definiowanego jako głosy nieważne / (głosy ważne + nieważne)) w wyborach samorządowych w 2014. Pierwszy histogram dotyczy całej Polski (27455 komisji), drugi województwa pomorskiego (1856) a trzeci Mazowieckiego (3574).

#!/usr/bin/Rscript
# Skrypt wykreśla histogramy dla danych z pliku ws2014_komisje.csv
# (więcej: https://github.com/hrpunio/Data/tree/master/ws2014_pobranie_2018)
#
par(ps=6,cex=1,cex.axis=1,cex.lab=1,cex.main=1.2)
komisje <- read.csv("ws2014_komisje.csv", sep = ';',
       header=T, na.string="NA");

komisje$ogn <- komisje$glosyNiewazne  / (komisje$glosy + komisje$glosyNiewazne) * 100;

summary(komisje$glosyNiewazne); fivenum(komisje$glosyNiewazne);
sX <- summary(komisje$ogn);
sF <- fivenum(komisje$ogn);
sV <- sd(komisje$ogn, na.rm=TRUE)
skewness <- 3 * (sX[["Mean"]] - sX[["Median"]])/sV

summary_label <- sprintf ("Śr = %.1f\nMe = %.1f\nq1 = %.1f\nq3 = %.1f\nW = %.2f", 
  sX[["Mean"]], sX[["Median"]], sX[["1st Qu."]], sX[["3rd Qu."]], skewness)

## ##
kpN <- seq(0, 100, by=2);
kpX <- c(0, 10,20,30,40,50,60,70,80,90, 100);
nn <- nrow(komisje)

h <- hist(komisje$ogn, breaks=kpN, freq=TRUE,
   col="orange", main=sprintf ("Rozkład odsetka głosów nieważnych\nPolska ogółem %i komisji", nn), 
   ylab="%", xlab="% nieważne", labels=F, xaxt='n' )
   axis(side=1, at=kpN, cex.axis=2, cex.lab=2)
   posX <- .5 * max(h$counts)
text(80, posX, summary_label, cex=1.4, adj=c(0,1))

## ##
komisje$woj <- substr(komisje$teryt, start=1, stop=2)

komisjeW <- subset (komisje, woj == "22"); ## pomorskie
nn <- nrow(komisjeW)
sX <- summary(komisjeW$ogn); sF <- fivenum(komisjeW$ogn);
sV <- sd(komisjeW$ogn, na.rm=TRUE)
skewness <- 3 * (sX[["Mean"]] - sX[["Median"]])/sV

summary_label <- sprintf ("Śr = %.1f\nMe = %.1f\nq1 = %.1f\nq3 = %.1f\nW = %.2f", 
  sX[["Mean"]], sX[["Median"]], sX[["1st Qu."]], sX[["3rd Qu."]], skewness)

h <- hist(komisjeW$ogn, breaks=kpN, freq=TRUE,
   col="orange", main=sprintf("Rozkład odsetka głosów nieważnych\nPomorskie %i komisji", nn), 
   ylab="%", xlab="% nieważne", labels=T, xaxt='n' )
   axis(side=1, at=kpX, cex.axis=2, cex.lab=2)
   posX <- .5 * max(h$counts)
text(80, posX, summary_label, cex=1.4, adj=c(0,1))

komisjeW <- subset (komisje, woj == "14"); ## mazowieckie
nn <- nrow(komisjeW)
sX <- summary(komisjeW$ogn); sF <- fivenum(komisjeW$ogn);
sV <- sd(komisjeW$ogn, na.rm=TRUE)
skewness <- 3 * (sX[["Mean"]] - sX[["Median"]])/sV

summary_label <- sprintf ("Śr = %.1f\nMe = %.1f\nq1 = %.1f\nq3 = %.1f\nW = %.2f", 
  sX[["Mean"]], sX[["Median"]], sX[["1st Qu."]], sX[["3rd Qu."]], skewness)

h <- hist(komisjeW$ogn, breaks=kpN, freq=TRUE,
   col="orange", main=sprintf("Rozkład odsetka głosów nieważnych\nMazowieckie %i komisji", nn), 
   ylab="%", xlab="% nieważne", labels=T, xaxt='n' )
   axis(side=1, at=kpX, cex.axis=2, cex.lab=2)
   posX <- .5 * max(h$counts)
text(80, posX, summary_label, cex=1.4, adj=c(0,1))

Wyniki są takie oto (indywidualne wykresy tutaj: #01 #02 #03):

Rozkłady odsetka poparcia dla PSL/PiS/PO w wyborach samorządowych w 2014 w całej Polsce, w miastach/poza miastami oraz w poszczególnych województwach. Poniższy skrypt generuje łącznie 60 wykresów słupkowych:

#!/usr/bin/Rscript
# Skrypt wykreślna różnego rodzaju histogramy dla danych z pliku ws2014_komitety_by_komisja_T.csv
# (więcej: https://github.com/hrpunio/Data/tree/master/ws2014_pobranie_2018)
#
showVotes <- function(df, x, co, region, N, minN) {
   ## showVotes = wykreśla histogram dla województwa (region)
   kN <- nrow(df)
   sX <- summary(df[[x]], na.rm=TRUE);
   sV <- sd(df[[x]], na.rm=TRUE)
   ## współczynnik skośności Pearsona
   skewness <- 3 * (sX[["Mean"]] - sX[["Median"]])/sV

   summary_label <- sprintf ("Śr = %.1f\nMe = %.1f\nq1 = %.1f\nq3 = %.1f\nS = %.1f\nW = %.2f", 
     sX[["Mean"]], sX[["Median"]],
     sX[["1st Qu."]], sX[["3rd Qu."]], sV, skewness)

   if (minN < 1) {
   t <- sprintf("Rozkład głosów na %s\n%s ogółem %d komisji", co, region, kN ) } 
   else { t <- sprintf("Rozkład głosów za %s\n%s ogółem %d komisji (N>%d)", co, region, kN, minN ) } 

   h <- hist(df[[x]], breaks=kpN, freq=TRUE, col="orange", main=t, 
     ylab="%", xlab="% poparcia", labels=F, xaxt='n' )
     axis(side=1, at=kpN, cex.axis=2, cex.lab=2)
   ## pozycja tekstu zawierającego statystyki opisowe
   posX <- .5 * max(h$counts)
   text(80, posX, summary_label, cex=1.4, adj=c(0,1))
}

## Wczytanie danych; obliczenie podst. statystyk:
komisje <- read.csv("ws2014_komitety_by_komisja_T.csv", 
   sep = ';', header=T, na.string="NA");

komisje$ogn <- komisje$glosyNiewazne  / (komisje$glosy 
   + komisje$glosyNiewazne) * 100;

summary(komisje$PSL); summary(komisje$PiS); summary(komisje$PO);
fivenum(komisje$PSLp); fivenum(komisje$PiSp); fivenum(komisje$POp);

## ## ###
par(ps=6,cex=1,cex.axis=1,cex.lab=1,cex.main=1.2)
kpN <- seq(0, 100, by=2);
kpX <- c(0, 10,20,30,40,50,60,70,80,90, 100);
kN <- nrow(komisje)
region <- "Polska"
minTurnout <- 0

## cała Polska:
showVotes(komisje, "PSLp", "PSL", region, kN, minTurnout);
showVotes(komisje, "PiSp", "PiS", region, kN, minTurnout);
showVotes(komisje, "POp",  "PO",  region, kN, minTurnout);

## Cała Polska (bez małych komisji):
## ( późniejszych analizach pomijane są małe komisje)
minTurnout <- 49
komisje <- subset (komisje, glosyLK > minTurnout); 
kN <- nrow(komisje)

showVotes(komisje, "PSLp", "PSL", region, kN, minTurnout);
showVotes(komisje, "PiSp", "PiS", region, kN, minTurnout);
showVotes(komisje, "POp",  "PO",  region, kN, minTurnout);

## Typ gminy U/R (U=gmina miejska ; R=inna niż miejska)
komisjeW <- subset (komisje, typ == "U"); 
kN <- nrow(komisjeW)
region <- "Polska/g.miejskie"
showVotes(komisjeW, "PSLp", "PSL", region, kN, minTurnout);
showVotes(komisjeW, "PiSp", "PiS", region, kN, minTurnout);
showVotes(komisjeW, "POp",  "PO",  region, kN, minTurnout);

komisjeW <- subset (komisje, typ == "R"); 
kN <- nrow(komisjeW)
region <- "Polska/g.niemiejskie"
showVotes(komisjeW, "PSLp", "PSL", region, kN, minTurnout);
showVotes(komisjeW, "PiSp", "PiS", region, kN, minTurnout);
showVotes(komisjeW, "POp",  "PO",  region, kN, minTurnout);

## woj = dwucyfrowy kod teryt województwa:
komisje$woj <- substr(komisje$teryt, start=1, stop=2)

cN <- c("dolnośląskie", "dolnośląskie", "kujawsko-pomorskie",
 "lubelskie", "lubuskie", "łódzkie", "małopolskie", "mazowieckie",
 "opolskie", "podkarpackie", "podlaskie", "pomorskie", "śląskie",
 "świętokrzyskie", "warmińsko-mazurskie", "wielkopolskie",
 "zachodniopomorskie");
cW <- c("02", "04", "06", "08", "10", "12", "14", "16", "18",
 "20", "22", "24", "26", "28", "30", "32");

## wszystkie województwa po kolei:
for (w in 1:16) {
  wojS <- cW[w]
  ###region <- cN[w];
  region <- sprintf ("%s (%s)", cN[w], wojS);

  komisjeW <- subset (komisje, woj == wojS); ##

  showVotes(komisjeW, "PSLp", "PSL", region, kN, minTurnout);
  showVotes(komisjeW, "PiSp", "PiS", region, kN, minTurnout);
  showVotes(komisjeW, "POp",  "PO",  region, kN, minTurnout);
}
## ## koniec

Dla całej Polski wyniki są następujące:

Indywidualne wykresy zaś tutaj: #01 #02 #03 #04 #05 #06 #07 #08 #09 #10 #11 #12 #13 #14 #15 #16 #17 #18 #19 #20 #21 #22 #23 #24 #25 #26 #27 #28 #29 #30 #31 #32 #33 #34 #35 #36 #37 #38 #39 #40 #41 #42 #43 #44 #45 #46 #47 #48 #49 #50 #51 #52 #53 #54 #55 #56 #57 #58 #59 #60):

Się tak rozochociłem, że zapisałem się na Żuławy w Koło 2018 (99 PLN). W ramach przygotowań, że tak powiem taktycznych postanowiłem rozpoznać możliwości przeciwnika:-) Konkretnie ustalić jak jechali ci co się zapisali, a co już startowali w ŻwK w roku 2017 albo 2016. Zadanie zatem polega na odszukaniu na liście zgłoszeń tych co się zapisali na edycję 2018 i jednocześnie ukończyli ŻwK w latach 2016/2017. Oczywiście nie ręcznie, tylko automatem:

## poniższe ściąga plik z listą zapisanych
wget 'http://www.czasomierzyk.pl/zapisy2016/zulawywkolo/index.php?akcja=lista' -O ZwK2018.out

Plik HTML ma tak prostą strukturę, że jego zamiana (za pomocą wyrażeń regularnych) na CSV jest banalna. Jak już mam ten plik CSV, to porównuję go do połączonych wyników z lat 2017/2016 (też w formacie CSV). Skrypt mam co porównuje pliki CSV:

perl join_csvs.pl -fn1 ZwK201809190908.csv  -fs1 1,2 -fn2 ZwK16_17.csv -fs2 1,2

Porównuje pliki ZwK201809190908.csv oraz ZwK16_17.csv, w oparciu o (wspólną) wartości dla kolumn nr 1 oraz nr 2 (w tym przypadku są to kolumny zawierające nazwisko i imię). Innymi słowy fs1 c1,c2..., to klucz główny, a fs2 c1,c2, to klucz obcy. Skrypt wypisuje połączone wiersze odpowiadające tym wierszom dla, których klucz główny = klucz obcy. Na dziś (19 września) takich wierszy wypisał 55, (na 104 zgłoszenia na dystansie 140km), ale pomijam tych co startowali kiedyś na najkrótszym dystansie lub tych, którzy startowali wprawdzie na najkrótszym, ale mieli średnią mniejszą niż 24kmh (odpada w ten sposób 10 zostaje 45). Na koniec plik jest zapodawany do prostego skryptu rysującego wykres słupkowy:

z <- read.csv("ZwK_2018_vs_2017.csv", sep = ';',
   header=T, na.string="NA", dec=".");

s140 <- summary(z$speed)

z <- subset (z, ( speed > 16.0 )); ## bez maruderów

# wykres słupkowy
h <- hist(z$speed,
  breaks=c(18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35),
  freq=TRUE,
  col="orange",
  main="Dystans: 140 (biorący udział w latach 2017-16)",
  xlab="Prędkość średnia w latach 2017--16 [kmh]",ylab="L.kolarzy",
  labels=T, xaxt='n' )
  axis(side=1, at=c(18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35))
  text(38, 37, summary_label, cex = .8, adj=c(1,1) )

Jak widać paru ludków w okolicach 30kmh jest. Będzie za kim jechać.

Ze śladu GPX prostym skryptem wyciągam co trzeba tworząc plik CSV o następującej zawartości (nazwy kolumn: data-czas,wysokość,prędkość,dystans przebyty):

daytime;ele;speed;dist  

Teraz poniższym skryptem rysuję profile wysokości (wysokość/prędkość vs czas oraz wysokość/prędkość vs dystans)

library(reshape)
require(ggplot2)

graphWd <- 6
graphHt <- 5

args = commandArgs(trailingOnly = TRUE);

if (length(args)==0) { stop("Podaj nazwę pliku CSV", call.=FALSE) }

fileBase <- gsub(".csv", "", args[1]);
outFile1 <- paste (fileBase, "_1.pdf", sep = "");
outFile2 <- paste (fileBase, "_2.pdf", sep = "");

what <- args[2];

# http://stackoverflow.com/questions/7381455/filtering-a-data-frame-by-values-in-a-column
d <- read.csv(args[1], sep = ';',  header=T, na.string="NA");
coeff <- median(d$ele)/median(d$speed)
d$speed <- d$speed * coeff


p1 <- ggplot(d, aes(x = as.POSIXct(daytime, format="%Y-%m-%dT%H:%M:%SZ"))) +
  geom_line(aes(y = ele, colour = 'wysokość', group = 1), size=1.5) +
  geom_line(aes(y = speed, colour = 'prędkość', group = 1), size=.5) +
  stat_smooth(aes(y=speed, x=as.POSIXct(daytime, format="%Y-%m-%dT%H:%M:%SZ"), colour ='prędkość wygładzona')) +
  ylab(label="Wysokość [mnpm]") +
  xlab(label="czas") +
  scale_y_continuous( sec.axis = sec_axis(name="Prędkość [kmh]",  ~./ coeff)) +
  labs(colour = paste( what )) +
  theme(legend.position="top") +
  theme(legend.text=element_text(size=12));
p1
ggsave(file=outFile1, width=graphWd, height=graphHt )

p2 <- ggplot(d, aes(x = dist)) +
  geom_line(aes(y = ele, colour = 'wysokość', group = 1), size=1.5) +
  geom_line(aes(y = speed, colour = 'prędkość', group = 1), size=.5) +
  ##geom_smooth() +
  stat_smooth(aes(y=speed, x=dist, colour ='prędkość wygładzona')) +
  ylab(label="Wysokość [mnpm]") +
  xlab(label="dystans") +
  scale_y_continuous( sec.axis = sec_axis(name="Prędkość [kmh]",  ~./ coeff)) +
  labs(colour = paste( what )) +
  theme(legend.position="top") +
  theme(legend.text=element_text(size=12));
p2

ps <- stat_smooth(aes(y=speed, x=dist));

ggsave(file=outFile2, width=graphWd, height=graphHt )

Teraz na koniec ciekawostka. Mój smartfon produkuje pliki GPX z superdokładnym stemplem czasu np. 2018-08-23T04:52:43.168Z, na czym wysypuje się R. Po prostu usuwam część po kropce dziesiętnej oraz samą kropkę (tj. .168Z) i działa.