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