Weblog Tomasza Przechlewskiego [Zdjęcie T. Przechlewskiego]


scrum
random image [Photo gallery]
Zestawienie tagów
1-wire | 18b20 | 1wire | 2140 | 3rz | adamowicz | alsamixer | amazon | 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 | budyń | bursztyn | canon | cedewu | chaos | chello | chiller | chillerpl | chown | chujowetaśmy | ciasto | cmentarz | contour | coronavirus | covi19 | covid19 | cron | css | csv | cukinia | curl | cycling | d54250wykh | dbi | debian | dejavu | dhcp | dht22 | dia | docbook | dom | dp1500 | ds18b20 | duda | dulkiewicz | 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 | foto | france | francja | fripp | froggit | fuczki | fuji | fuse | gammu | garden | garmin | gawk | gazwyb | gdańsk | gdynia | gender | geo | geocoding | georgia | gft | ggplot | git | github | gmail | gmaps | gnokii | gnus | google | googlecl | googleearth | googlemaps | gotowanie | gphoto | gphoto2 | gps | gpsbabel | gpsphoto | gpx | gpx-viewer | greasemonkey | gruzja | grzyby | gw1000 | haldaemon | handbrake | hhi | historia | history | hitler | holocaust | holokaust | hp1000se | hpmini | humour | iblue747 | ical | iiyama | ikea | imagemagick | imap | inkscape | inne | internet | j10i2 | javascript | jhead | k800i | kajak | kamera | karob | kleinertest | kml | kmobiletools | knuth | kociewie kołem | kod | kolibki | komorowski | konwersja | krutynia | krynki | kuchnia | kurski | latex | latex2rtf | latex3 | lcd | legend | lenny | lesund | lewactwo | lgbt-folly | liberation | linksys | linux | lisp | lisrel | litwa | lizbona | logika | ltr | lubowla | lwp | lwów | m2wś | malta | mapquest | mapsource | marchew | marimekko | marvell | math | mathjax | mazury | mbank | mediolan | mencoder | mevo | mh17 | michalak | michlmayr | microsoft | monitor | mp4box | mplayer | ms | msc | mssql | msw | mswindows | mtkbabel | museum | muzyka | mymaps | mysql | nafisa | nanopi | natbib | navin | 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 | overclocking | ozbekiston | panoramio | paryż | pdf | pdfpages | pdftex | pdftk | pedophilia | perl | photo | photography | picasa | picasaweb | pim | pine | pis | pit | plotly | pls | plugin | po | 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 | rpi | rsync | rtf | ruby | rugby | rumunia | russia | rwc | rwc2007 | rwc2011 | rwc2019 | rzym | samba | sds011 | selenium | sem | sernik | sheevaplug | sienkiewicz | signature | 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 | wh2080 | wiedeń | wikicommons | wilno | win10 | windows | windows8 | wine | wioślarstwo | 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 | zakopane | zakupy | zawodzie | zdf | zdrowie | łeba | świdnica | żywność
Archiwum
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
Zmarli/zakażeni i zmarli/1mln




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.

url | Thu, 24/09/2020 05:35 | tagi: , , ,
Jeszcze jeden skrypt do wizualizacji danych nt COVID19

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?

url | Thu, 17/09/2020 20:37 | tagi: , , ,
Rosyjskie dane dotyczące urodzin i zgonów






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 )
url | Wed, 16/09/2020 09:37 | tagi: , , ,
Wydłubywanie danych nt/ COVID19 z tweetów MZ

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

url | Wed, 02/09/2020 21:46 | tagi: , , , ,
Worldometer vs ECDC

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:





url | Sun, 26/04/2020 05:17 | tagi: , ,
Wględne tempo wzrostu (koronowirusa)

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

url | Tue, 31/03/2020 05:16 | 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: , , ,
Dzienne dane dot. wypadków drogowych w Polsce

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

url | Wed, 25/03/2020 07:29 | tagi: , ,
Mevo: koncentracja rowerów na stacjach

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

url | Thu, 01/08/2019 12:38 | tagi: , , , ,
Wybory 2018. Różnica w liczbie mandatów do Sejmików

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ą).

url | Fri, 26/10/2018 05:47 | tagi: , , ,
Wybory wójtów/burmistrzów/prezydentów

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

url | Thu, 18/10/2018 17:31 | tagi: , , ,
Koniec pobierania danych wyborczych

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

url | Thu, 11/10/2018 08:31 | tagi: , , ,
Wybory 2014 i jeszcze więcej rozkładów

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

url | Tue, 02/10/2018 17:08 | tagi: , , ,
Żuławy w Koło 2018

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

url | Wed, 19/09/2018 12:54 | tagi: , , ,
Rysowanie profilu wysokości w R

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.

url | Fri, 31/08/2018 07:40 | tagi: , ,
time plot tygodniowej liczby twitów

Załóżmy, że plik CSV zawiera liczbę opublikowanych twitów (dane tygodniowe). Problem: przedstawić szereg w postaci przebiegu czasowego (time plot). Taki skrypt R wymyśliłem do zrealizowania tego zadania:

require(ggplot2)

args <- commandArgs(TRUE)
ttname <- args[1];
file <- paste(ttname, ".csv", sep="")
filePDF <- paste(ttname, ".pdf", sep="")

d <- read.csv(file, sep = ';',  header=T, na.string="NA", );
## Plik CSV jest postaci:

##str(d)

## wiersze 1,2 + ostatni są nietypowe (usuwamy)
rows2remove <- c(1, 2, nrow(d));
d <- d[ -rows2remove, ];

## szacujemy prosty model trendu
lm <- lm(data=d, posts ~ no ); summary(lm)
posts.stats <- fivenum(d$posts);
posts.mean <- mean(d$posts);
sumCs <- summary(d$posts);

otherc <- coef(lm);
# W tytule średnia/mediana i równanie trendu
title <- sprintf ("Weekly for %s # me/av = %.1f/%.1f (y = %.2f x + %.1f)", 
  ttname, sumCs["Median"], sumCs["Mean"], otherc[2], otherc[1] );

##str(d$no)
## Oś x-ów jest czasowa
## Skróć yyyy-mm-dd do yy/mmdd
d$date <- sub("-", "/", d$date) ## zmienia tylko pierwszy rr-mm-dd
d$date <- sub("-", "", d$date) ## usuwa mm-dd
d$date <- gsub("^20", "", d$date) ## usuwa 20 z numeru roku 2018 -> 18
weeks <- length(d$no);
## https://stackoverflow.com/questions/5237557/extract-every-nth-element-of-a-vector
## Na skali pokaż do 20 element /dodaj ostatni `na pałę' (najwyżej zajdą na siebie)
## możnaby to zrobić bardziej inteligentnie ale nie mam czasu
scaleBreaks <- d$no[c(seq(1, weeks, 20), weeks)];
scaleLabs <- d$date[c(seq(1, weeks, 20), weeks)];

ggplot(d, aes(x = no, y = posts)) +
  geom_line() +
  ggtitle(title) +
  ylab(label="#") +
  xlab(label=sprintf("time (yy/mmdd) n=%d", weeks )) +
  scale_x_continuous(breaks=scaleBreaks, labels=scaleLabs) +
  geom_smooth(method = "lm")

ggsave(file=filePDF)  

url | Wed, 21/02/2018 13:56 | tagi: , ,
Żuławy wKoło 2017 -- podsumowanie



Podsumowanie wyników dla lat 2015--2017

z16 <- read.csv("wyniki_zulawy_2016_D.csv", sep = ';',  header=T, na.string="NA", dec=",");
aggregate (z16$time, list(Numer = z16$dist), summary)
z16$year <- 2016;

z15 <- read.csv("wyniki_zulawy_2015_D.csv", sep = ';',  header=T, na.string="NA", dec=",");
aggregate (z15$time, list(Numer = z15$dist), summary)
z15$year <- 2015;

z17 <- read.csv("wyniki_zulawy_2017_D.csv", sep = ';',  header=T, na.string="NA", dec=".");
aggregate (z17$time, list(Numer = z17$dist), summary)
z17$year <- 2017;

zz15 <- z15[, c("dist", "kmH", "time", "year")];
zz16 <- z16[, c("dist", "kmH", "time", "year")];
zz17  <- z17[, c("dist", "kmH", "time", "year")];

zz <- rbind (zz15, zz16, zz17);

## tylko dystans 140
zz140 <- subset (zz, ( dist == 140 ));
sum140 <- aggregate (zz140$kmH, list(Numer = zz140$year), summary)

boxplot (kmH ~ year, zz140, ylab = "Śr.prędkość [kmh]", col = "yellow", main="140km" )

## tylko dystans 55
zz75 <- subset (zz, ( dist > 60 & dist < 90 ));
sum75 <- aggregate (zz75$kmH, list(Numer = zz75$year), summary)
sum75
boxplot (kmH ~ year, zz75, ylab = "Śr.prędkość [kmh]", col = "yellow", main="80/75km" )

## tylko dystans 55
zz55 <- subset (zz, ( dist < 60 ));
sum55 <- aggregate (zz55$kmH, list(Numer = zz55$year), summary)
sum55
xl <- paste ("średnie 2015=", sum55$x[1,4], "kmh   2016=",
  sum55$x[2,4], "kmh   2017=", sum55$x[3,4], " kmh")

  boxplot (kmH ~ year, zz55, xlab = xl,
  ylab = "Śr.prędkość [kmh]", col = "yellow", main="55km" )

A ja (numer 418) byłem 70 w kategorii 140 km, z czasem 5:42:03 co dało 24,56 kmh przeciętną. Do pierwszego bufetu się spinałem, potem już nie...

url | Thu, 28/09/2017 05:17 | tagi: , ,
Chancellor Merkel victory for a visual person

Change in number of seats won by party (AfD is brown of course regardless official party colors :-):-)

library(ggplot2)

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

ggplot(df, aes(x=party, y=diff, fill=party )) +
geom_bar(stat="identity") +
geom_text(aes(label=diff), vjust=-0.5) +
labs(x = "", y="change") +

ggtitle("German elections results (#MP change)") +

## AfD is brown regardless official party colors :-)
scale_fill_manual(values=c("#8B4513", "#56B4E9",
"yellow", "green", "red", "#ff6666") )
url | Mon, 25/09/2017 08:43 | tagi: , ,
Przed Żuławy w Koło 2017

Jutro planuję przejechać 140km biorąc udział w imprezie pn. Żuławy wKoło 2017. Niby Żuławy a profil trasy sugeruje jakieś istotne wzniesienie w okolicach 25--40km:

Uważnie przyjrzenie się liczbom (zwłaszcza na osi OY) pozwala stwierdzić, że jest to złudzenie, wynikające z różnicy w jednostkach miary obu osi (kilometry vs metry). W rzeczywistości góra tam jest symboliczna o czym można się przekonać robiąc wykres nachyleń. Żeby pozbyć się przypadkowych błędów związanych z niedokładnością pomiaru oryginalne 673 punktowe dane zostały zmienione na 111 punktowe (uśrednienie minimum 1 km) lub 62 punktowe (uśrednienie minimum 2 km).

Przy czym uśrednienie minimum $x$ oznacza obliczenie nachylenia dla najkrótszego odcinka kolejnych $n$ punktów z oryginalnego śladu GPX, który będzie dłuższy niż $x$.

Skrypty R/dane są tutaj. Oryginale ślady GPX/TCX skopiowane ze strony ŻwK są tutaj.

url | Sat, 23/09/2017 17:02 | tagi: , , ,
Czytelnictwo prasy

Punktem wyjścia są dane ze strony ZKDP (w formacie Excel.) Ponieważ pobieram je od pewnego czasu mam tego więcej niż jest na ww. stronie bo od stycznia 2015 do lipca 2017, czyli 31 plików. Ręczna konwersja byłaby zatem ciut za bardzo czasochłonna.

  for i in *.xls do
    oocalc --headless --convert-to csv $i ;
    # albo ssconvert -v $i `basename $i .xls`.csv ;
    done
  # Wyciągam dane dotyczące sprzedaży ogółem dla SE
  grep 'Super Ex' *.csv | awk -F ',' '{print $7} ' > se_sales.csv
  # Analogicznie dla innych tytułów

Uwaga: program ssconvert znajduje się w pakiecie gnumeric, oocalc to oczywiście składni Libre/OpenOffice.

Wielkości sprzedaży dla trzech najpoczytniejszych tytułów pokazują wykresy liniowe (pierwszy w tys egz. a drugi w procentach nakładu ze stycznia 2015 r.)


Sprzedaż w tys egz.

Sprzedaż w % poziomu ze stycznia 2015
library(ggplot2)
library(reshape2)

df <- read.csv("newspaper_sales_2015-17.csv", sep = ';',
               header=T, na.string="NA");

meltdf <- melt(df,id="month")

ggplot(meltdf,aes(x=month, y=value, colour=variable, group=variable)) +
  geom_line() +
  ylab(label="sales [ths]") +
  theme(legend.title=element_blank()) +
  scale_x_discrete (breaks=c("2015-01-01", "2015-06-01",
     "2016-01-01", "2016-06-01", "2017-01-01",  "2017-06-01"),
  labels=c("2015-01", "2015-06", "2016-01", "2016-06",
     "2017-01", "2017-06")  )

# https://stackoverflow.com/questions/10085806/extracting-specific-columns-from-a-data-frame
obs <- df[,c("month")]

normalize <- function(x) { return (x /x[1] * 100 )  }
dfN <- as.data.frame(lapply(df[-1], normalize))

# https://stackoverflow.com/questions/10150579/adding-a-column-to-a-data-frame
dfN["month"] <- obs

str(dfN)

meltdf <- melt(dfN,id="month")

# https://www.r-bloggers.com/what-is-a-linear-trend-by-the-way/
pN <- ggplot(meltdf,
 aes(x=month, y=value, colour=variable, group=variable)) + geom_line() +
 ylab(label="sales [ths]") +
 theme(legend.title=element_blank()) +
 stat_smooth(method = "lm", se=F) +
  scale_x_discrete (breaks=c("2015-01-01", "2015-06-01",
     "2016-01-01", "2016-06-01", "2017-01-01",  "2017-06-01"),
  labels=c("2015-01", "2015-06", "2016-01",
  "2016-06", "2017-01", "2017-06")  )

  pN

Spadek widoczny na wykresach można określić liczbowo na przykład szacując linię trendu:

# Trend liniowy
# http://t-redactyl.io/blog/2016/05/creating-plots-in-r-using-ggplot2-part-11-linear-regression-plots.html

# http://r-statistics.co/Time-Series-Analysis-With-R.html
seq = c (1:nrow(dfN))
dfN["trend"] <- seq

trendL.gw <- lm(data=dfN, gw ~ trend )
trendL.fakt <- lm(data=dfN, fakt ~ trend )
trendL.se <- lm(data=dfN, se ~ trend )

trendL.gw
trendL.fakt
trendL.se

Współczynniki trendu dla GW, Faktu i SE są odpowiednio równe -1.114 -0.6415 -0.4301, co należy interpretować następująco: przeciętnie z miesiąca na miesiąc nakład spada o 1,11%, 0,64% oraz 0,43% nakładu ze stycznia 2015 r., odpowiednio dla GW, Faktu i SuperExpresu.

Dane i wyniki są tutaj

url | Mon, 11/09/2017 17:33 | tagi: , , ,
Żuławy w koło 2016

Żuławy w koło to maraton rowerowy (czyli przejazd rowerem na dłuższym dystansie -- nie mylić z wyścigiem) organizowany od paru lat na Żuławach jak nazwa wskazuje. Sprawdziłem jak ta impreza wyglądała pod kątem prędkości w roku 2016. W tym celu ze strony Wyniki żUŁAWY wKOŁO 2016 ściągnąłem stosowny plik PDF z danymi, który następnie skonwertowałem do pliku w formacie XLS (Excel) wykorzystując konwerter on-line tajemniczej firmy convertio.pl. Tajemniczej w tym sensie, że nie znalazłem informacji kto i po co tą usługę świadczy.

Konwersja (do formatu CSV) -- jak to zwykle konwersja -- nie poszła na 100% poprawnie i wymagała jeszcze circa 30 minutowej ręcznej obróbki. Być może zresztą są lepsze konwertery, ale problem był z gatunku banalnych i wolałem stracić 30 minut na poprawianiu wyników konwersji niż 2 godziny na ustalaniu, który z konwerterów on-line konwertuje ten konkretny plik PDF (w miarę) bezbłędnie.

Po konwersji wypadało by sprawdzić (chociaż zgrubnie) czy wszystko jest OK.

## Czy każdy wiersz zawieraja 9 pól (powinien)
$ awk -F ';' 'NF != 9 {print NR, NF}' wyniki_zulawy_2016S.csv

## Ilu było uczestników na dystansie 140km?
$ awk -F ';' '$7 ==140 {print $0}' wyniki_zulawy_2016S.csv | wc -l
133

## Ilu było wszystkich (winno być 567 + 1 nagłówek)
$ cat wyniki_zulawy_2016S.csv | wc -l
568 # ok!

Przykładowy wykres pudełkowy

Do analizy statystycznej wykorzystano wykres pudełkowy (porównanie wyników na różnych dystansach) oraz histogram (rozkład średnich prędkości na dystansie 140km). BTW gdyby ktoś nie wiedział co to jest wykres pudełkowy to wyjaśnienie jest na rysunku obok. Objaśnienie: Me, $Q_1$, $Q_3$ to odpowiednio mediana i kwartyle. Dolna/górna krawędź prostokąta wyznacza zatem rozstęp kwartylny (IQR). Wąsy ($W_L$/$W_U$) są wyznaczane jako 150% wartości rozstępu kwartylnego. Wartości leżące poza ,,wąsami'' (nietypowe) są oznaczane kółkami.

Ww. wykresy wygenerowano następującym skryptem:

#
co <- "Żuławy wKoło 2016"
#
z <- read.csv("wyniki_zulawy_2016_C.csv", sep = ';',
  header=T, na.string="NA", dec=",");

aggregate (z$meanv, list(Numer = z$dist), fivenum)

boxplot (meanv ~ dist, z, xlab = "Dystans [km]",
    ylab = "Śr.prędkość [kmh]", col = "yellow", main=co )

## tylko dystans 140
z140 <- subset (z, ( dist == 140 ));

## statystyki zbiorcze
s140 <- summary(z140$meanv)
names(s140)

summary_label <- paste (sep='', "Średnia = ", s140[["Mean"]], 
  "\nMediana = ", s140[["Median"]],
  "\nQ1 = ", s140[["1st Qu."]],  "\nQ3 = ", s140[["3rd Qu."]],
  "\n\nMax = ", s140[["Max."]] )
# drukuje wartości kolumny meanv
# z140$meanv
# drukuje wartości statystyk zbiorczych
s140

# wykres słupkowy
h <- hist(z140$meanv, breaks=c(14,18,22,26,30,34,38), freq=TRUE, 
   col="orange", 
   main=paste (co, "[140km]"), # tytuł
   xlab="Prędkość [kmh]",ylab="L.kolarzy", labels=T, xaxt='n' )
# xaxt usuwa domyślną oś 
# axis definiuje lepiej oś OX
axis(side=1, at=c(14,18,22,26,30,34,38))
text(38, 37, summary_label, cex = .8, adj=c(1,1) )

Dane i wyniki są tutaj

url | Mon, 11/09/2017 08:10 | tagi: , , ,