Weblog Tomasza Przechlewskiego [Zdjęcie T. Przechlewskiego]


scrum
random image [Photo gallery]
Zestawienie tagów
1-wire | 18b20 | 1wire | 2140 | 3rz | adamowicz | afera | alsamixer | amazon | amber | amman | anniversary | antypis | apache | api | applebaum | arm | armenia | astronomy | asus | atom.xml | awk | aws | bachotek | bakłażan | balcerowicz | balta | banan | bash | batumi | berlin | białowieża | białystok | bibtex | bieszczady | biznes | blogger | blogging | blosxom | bme280 | bono | borne-sulinowo | breugel | bt747 | budapeszt | budyniowo | budyń | bursztyn | campagnolo | canon | cedewu | chaos | chello | chiller | chillerpl | chown | christophe dominici | chujowetaśmy | ciasto | cmentarz | contour | coronavirus | covi19 | covid | covid19 | cron | css | csv | cukinia | curl | cycling | d54250wykh | darkages | dbi | debian | dejavu | dhcp | dht22 | dia | docbook | dom | dp1500 | ds18b20 | duda | dulkiewicz | dulkiewiczowa | dyndns | dynia | ebay | economy | ecowitt | ekonomia | elka | elm | emacs | emacs23 | english | ep | erasmus | erasmusplus | ess | eu | eurostat | excel | exif | exiftool | f11 | fc | fc11 | fc15 | fc29 | fc5 | fc8 | fedora | fedora21 | fenix | ffmpeg | finepix | firefox | flickr | folau | fontforge | fontspec | fonty | food | fop | forms | foto | france | francja | fripp | froggit | fuczki | fuji | fuse | gammu | garden | garmin | gas | gawk | gazwyb | gdańsk | gdynia | gender | geo | geocoding | georgia | gft | ggplot | ghost | git | github | gmail | gmaps | gnokii | gnus | google | google apps script | googlecl | googleearth | googlemaps | gotowanie | gphoto | gphoto2 | gps | gpsbabel | gpsphoto | gpx | gpx-viewer | greasemonkey | gruzja | grzyby | gus | gw1000 | haldaemon | handbrake | hhi | historia | history | hitler | holocaust | holokaust | hp1000se | hpmini | humour | iblue747 | ical | iiyama | ikea | imagemagick | imap | inkscape | inne | internet | j10i2 | javascript | jhead | jordania | k800i | kajak | kamera | karob | kibbeh | kleinertest | kml | kmobiletools | knuth | kociewie kołem | kod | kolibki | komorowski | konwersja | krutynia | krynki | kuchnia | kurski | kłamstwo | latex | latex2rtf | latex3 | lcd | legend | lenny | lesund | lewactwo | lgbt-folly | liban | liberation | linksys | linux | lisp | lisrel | litwa | lizbona | logika | ltr | lubowla | lwp | lwów | m2wś | malta | mapquest | mapsource | maradona | marchew | marimekko | marvell | math | mathjax | mazury | mbank | mediolan | mencoder | mevo | mex | mh17 | michalak | michlmayr | microsoft | monitor | mp4box | mplayer | ms | msc | mssql | msw | mswindows | mtkbabel | museum | muzyka | mymaps | mysql | mz | nafisa | nanopi | natbib | navin | neapol | nekrolog | neo | neopi | netbook | niemcy | niemieckie zbrodnie | nikon | nmea | nowazelandia | nuc | nxml | oauth | oauth2 | obituary | ocr | odessa | okular | olympus | ooffice | ooxml | opera | osm | otf | otftotfm | other | ov5647 | overclocking | ozbekiston | padwa | panoramio | paryż | pdf | pdfpages | pdftex | pdftk | pedophilia | perl | photo | photography | pi | picasa | picasaweb | pim | pine | pis | pit | pizero | plain | plotly | pls | plugin | po | podcast | podlasie | podróże | pogoda | politics | polityka | polsat | portugalia | postęp | powerpoint | połtawa | prelink | problem | propaganda | pseudointeligencja | pstoedit | putin | python | pywws | r | r1984 | radio | random | raspberry | raspberry pi | raspberrypi | raspbian | refugees | relaxng | ridley | router | rower | rowery | roztocze | rpi | rsync | rtf | ruby | rugby | rumunia | russia | rwc | rwc2007 | rwc2011 | rwc2019 | rzym | salerno | samba | sds011 | selenium | sem | senah | sernik | sheevaplug | sienkiewicz | signature | sikorski | sks | skype | skytraq | smoleńsk | sqlite | srtm | sshfs | ssl | staszek wawrykiewicz | statistcs | statistics | stats | statystyka | stix | stretch | supraśl | suwałki | svg | svn | swanetia | swornegacie | szwajcaria | słowacja | tbilisi | terrorism | tesseract | tex | texgyre | texlive | thunderbird | tomato | totalnaopozycja | tourism | tramp | trang | transylwania | truetype | trzaskowski | ttf | turcja | turkey | turystyka | tusk | tv | tv5monde | tweepy | twitter | tykocin | typetools | ubuntu | uchodźcy | udev | ue | ukraina | umap | unix | upc | updmap | ups | utf8 | uzbekistan | varia | video | vienna | virb edit | virbedit | vostro | wammu | wdc | wdfs | weather | weathercloud | webcam | webdav | webscrapping | weewx | wenecja | wh2080 | wiedeń | wikicommons | wilno | win10 | windows | windows8 | wine | wioślarstwo | wojna | word | wordpress | wrt54gl | ws1080 | wtyczka | wunderground | ww2 | www | wybory | wybory2015 | włochy | węgry | xemex | xetex | xft | xhtml | xine | xml | xmllint | xsd | xslt | xvidtune | youtube | yum | zaatar | zakopane | zakupy | zawodzie | zdf | zdrowie | zeropi | zgarden | zgony | zprojekt | łeba | świdnica | żywność
Archiwum
02/2023 | 01/2023 | 11/2022 | 10/2022 | 09/2022 | 07/2022 | 06/2022 | 04/2022 | 03/2022 | 02/2022 | 12/2021 | 09/2021 | 03/2021 | 01/2021 | 12/2020 | 11/2020 | 10/2020 | 09/2020 | 08/2020 | 07/2020 | 04/2020 | 03/2020 | 02/2020 | 01/2020 | 12/2019 | 11/2019 | 10/2019 | 09/2019 | 08/2019 | 07/2019 | 06/2019 | 04/2019 | 02/2019 | 01/2019 | 12/2018 | 11/2018 | 10/2018 | 09/2018 | 08/2018 | 07/2018 | 05/2018 | 04/2018 | 03/2018 | 02/2018 | 01/2018 | 11/2017 | 10/2017 | 09/2017 | 08/2017 | 07/2017 | 06/2017 | 05/2017 | 04/2017 | 03/2017 | 02/2017 | 01/2017 | 12/2016 | 11/2016 | 10/2016 | 09/2016 | 08/2016 | 06/2016 | 05/2016 | 04/2016 | 02/2016 | 12/2015 | 11/2015 | 09/2015 | 07/2015 | 06/2015 | 05/2015 | 02/2015 | 01/2015 | 12/2014 | 09/2014 | 07/2014 | 06/2014 | 04/2014 | 02/2014 | 01/2014 | 12/2013 | 11/2013 | 10/2013 | 09/2013 | 08/2013 | 07/2013 | 05/2013 | 04/2013 | 03/2013 | 02/2013 | 01/2013 | 12/2012 | 11/2012 | 10/2012 | 09/2012 | 08/2012 | 07/2012 | 05/2012 | 03/2012 | 02/2012 | 01/2012 | 12/2011 | 11/2011 | 10/2011 | 09/2011 | 08/2011 | 07/2011 | 06/2011 | 05/2011 | 04/2011 | 03/2011 | 02/2011 | 01/2011 | 12/2010 | 11/2010 | 10/2010 | 09/2010 | 08/2010 | 07/2010 | 06/2010 | 05/2010 | 04/2010 | 03/2010 | 02/2010 | 01/2010 | 12/2009 | 11/2009 | 10/2009 | 09/2009 | 08/2009 | 07/2009 | 06/2009 | 05/2009 | 04/2009 | 03/2009 | 02/2009 | 01/2009 | 12/2008 | 11/2008 | 10/2008 | 09/2008 | 08/2008 | 07/2008 | 06/2008 | 05/2008 | 04/2008 | 03/2008 | 02/2008 | 01/2008 | 12/2007 | 11/2007 | 10/2007 | 09/2007 | 08/2007 | 07/2007 |
O stronie
wykorzystywany jest blosxom plus następujące wtyczki: tagging, flatarchives, rss10, lastbuilddatexhtmlmime. Niektóre musiałem dopasować nieco do swoich potrzeb. Więcej o blosxom jest tutaj
Subskrypcja
RSS 1.0
Zgony nadmiarowe w Europie

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

url | Fri, 10/12/2021 05:56 | tagi: , , , ,
Dane Eurostatu nt zgonów/urodzeń

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)

url | Wed, 25/03/2020 08:51 | tagi: , , ,