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