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