Weblog Tomasza Przechlewskiego

Darrell Huff. How to lie with statistics

Bill Gates poleca w 2015

How to lie with statistics

Zdjęcie Gatesa (z 2015 roku) w połączeniu z faktem, że Gates finasował badania w dziedzinie epidemiologii (na John Hopkins University) stało się ,,dowodem'' dla różnych szurów, których w USA nie brakuje, iż za pandemią COVID19 stał Gates.

A book written by Darrell Huff in 1954 presenting an introduction to statistics for the general reader. Not a statistician, Huff was a journalist [...]

In the 1960/1970s, it became a standard textbook introduction to the subject of statistics for many college students [...] one of the best-selling statistics books in history.

https://en.wikipedia.org/wiki/How_to_Lie_with_Statistics

Książeczka składa się z 10 rozdziałów i jest napisana w prowokacyjny, sposób (nienaukowy). Nie była przetłumaczona na język polski. Poszczególne rozdziały można powiedzieć przeszły do legendy i jak się wpisze tytuł rozdziału do google to zwykle można znaleźć setki tysięcy stron cytujących... Osobiście nie widzę nic aż tak nadzwyczajnego w tej książce. Przedstawia kilkanaście sposobów manipulacji, w miarę oczywistych. Miejscami gubi wątek w tym sensie, że są rozdziały lepsze (zaznaczone plusem poniżej) i gorsze. Ale ponieważ jest tak znana to poniżej strzeszczenie:

r1+: a sample with the built-in bias czyli niereprezentatywność próby; że ciężko jest zebrać próbę reprezentatywną (z różnych powodów).

r2: the well chosen average czyli sztuczki nt. średniej. Zarówno co jest uśredniane (who's included), jak i jak jest uśredniane (średnia vs mediana)

r3+: the little figures that are not there. Niejasne/nieznane szczegóły wyników analizy (statystycznie nieistotne rezulataty ogłaszane bez podania, że są nieistotne--albo średnie dla rozkładów daleko różnych od normalnych)

r4: to samo co #r3 przy założeniu że pomiar jest mocno przybliżony przez co zaobserwowane różnice nie mają specjalnie znaczenia (bo ewentualny błąd jest większy niż różnice)

r5+: The gee-whiz graph aka zmyłkowe wykresy (głównie nie zaczynająca się od zera oś 0Y) (cf https://en.wikipedia.org/wiki/Gee_Whiz albo https://en.wikipedia.org/wiki/Misleading_graph)

r6+: The one dimensional picture aka zmyłkowe wykresy cd (porównywanie jednowymiarowych wielkości w 2D albo 3D; cf https://thejeshgn.com/2017/11/17/how-to-lie-with-graphs/)

r7+: semiattached figure. Using one thing as a way to claim proof of something else, even though there's no correlation between the two (teza i dowód nie są ze sobą powiązane niczym oprócz wrażenia że są; https://www.secjuice.com/the-semi-attached-figure/)

r8: post hoc Rides Again. Korelacja to nie przyczynowość; dla mnie najbardziej mętny rozdział ale też temat chyba najtrudniejszy do przybliżenia na poziomie Idiots Guide

r9: How to statisticulate: Misinforming people by the use of statistical material might be called statistical manipulation, in a word, Statisticulation. (ten rozdzialik to podsumowanie r1--r8)

r10++: how to talk back to statistics. Dwa plusy to nie przypadek bo chyba najciekawszy: Jak się nie dać oszukać kiepskiej statystyce w pięciu krokach.

Pięć kroków Huffa

Who Says So? (ludzie mają interesy, osoby zainteresowane mogą nie mówić prawdy);

How Does He Know? (pomiar jest często wysoce wadliwy);

What's Missing? (analiza jest niejasna/niepełna);

Many figures (liczb nie rysunków) lose meaning because a comparison is missing. Mój przykład: kobiety w PL nie rodzą dzieci; przeciętny wiek matki w momencie urodzenia dziecka to 27 lat. [czego NIE powiedziano: W całej Europie tak jest]

Did Somebody Change The Subject? (czy teza i dowód są logicznie powiązane czy tylko sprawiają takie wrażenie)

Does It Make Sense? (ogólnie czy coś z tego wynika na poziomie zdrowego rozsądku)

Bibliografia

Darrell Huff. How to lie with statistics (142 strony/a5) https://en.wikipedia.org/wiki/How_to_Lie_with_Statistics

url | Sun, 01/01/2023 14:57 | tagi: inne, statystyka

Zgony nadmiarowe w Europie

Że czwarta fala szaleje to nadmiarowe zgony w Europie sobie liczę, w której to kategorii Polska cytując poetę trzecie miejsce w świecie (Kazik Staszewski/Amnezja)

Dane pobieram za pomocą get_eurostat (pakiet eurostat):

library("eurostat")
# demo_r_mwk_ts = Deaths by week and sex
z <- get_eurostat("demo_r_mwk_ts",  stringsAsFactors = FALSE) %>%
  mutate (time = as.character(time)) %>%
  mutate (year = as.numeric(substr(time, 1, 4)), 
          week = as.numeric(substr(time, 6,7)),
          sex = as.factor(sex),
          geo = as.factor(geo)) %>%
  select (sex, geo, year, week, value=values) %>%
  filter (sex == 'T')

Końcowy filter usuwa liczbę zgonów dla kobiet/mężczyzn osobno, bo chcemy tylko analizować ogółem czyli Total.

W pliku ecdc_countries_names_codes3166.csv są nazwy krajów (bo w bazie Eurostatu są tylko ISO-kody):

## country names
nn <- read.csv("ecdc_countries_names_codes3166.csv", sep = ';',  header=T, na.string="NA" )

Teraz liczę 5-letnie średnie dla okresu przed pandemią czyli dla lat 2015--2019 włącznie:

## mean weekly deaths 2015--2019
z0 <- z %>% filter ( year >= 2015 & year < 2020) %>% 
  group_by(geo, week) %>% 
  summarise (d0 = mean(value, na.rm=T))

Oraz tygodniowe liczby zgonów dla lat 2020--2021

## weekly deaths 2020--2021
z1 <- z %>% filter ( year > 2019 ) %>% 
  group_by(geo, year, week) %>% 
  summarise ( d1 = sum(value))

Łączę z0 z z1; liczę różnice między zgonami 2020--2021 a średnimi. Najpierw usuwam wiersze z NA w kolumnach zawierających zgony/średnie (drop_na). To w szczegolności usunie wiersze z tygodni (w roku 2021), dla których jeszcze nie ma danych.

## join z0 z1 and compute differences
zz <- left_join(z0, z1, by=c("week", "geo")) %>% 
  drop_na(d0,d1) %>%
  mutate (exp = (d1 - d0)/d0 * 100, exm = d1 - d0 )
zz <- left_join(zz, nn, by=c('geo'='id'))

Osobno liczę numer ostatniego raportowanego tygodnia w roku 2021. Pytanie o numer ostatniego w ogóle jest bez sensu, bo w 2020 będzie to ostatni tydzień. Ten fragment ciut mało robust jest bo w 2022 trzeba go będzie zmodyfikować

## if NA then the country stop reporting in 2020
zz.last.week <- zz %>% filter (year == 2021) %>% 
  group_by (geo) %>% 
  summarise (lw = last(week)) %>%
  drop_na(lw)
  
## najbardziej aktualny raportowany tydzień (różne kraje raportują różnie)
latestweek <- max (zz.last.week$lw)

Uwaga: niektóre kraje nie raportują w 2021, więc ich nie będzie w ramce zz.last.week z uwagi na filter.

Ramka zzt zawiera sumy różnic; bezwzględne (exm) oraz względne w procentach względem poziomu średniego 2015--2019 (exp):

### total exmort 
zzt <- zz %>% group_by(geo) %>%
   summarise(country=last(country), 
   exm = sum(exm),
   nm = sum(d0)) %>%
  mutate (exmp = exm/nm * 100)

Dodajemy informację z ramki zz.last.week; jeżeli w tej ramce nie ma kraju (bo nie raportował w 2021) to usuwamy go (drop_na(lw)). Wypierniczamy kraje raportujące więcej niż 6 tygodni po kraju (krajach), który przysłały raport najpóźniej. Np jak jakiś kraj-prymus dostarczył w 47 tygodniu to tylko kraje raportujące z tygodni 41 i później zostają a inne są usuwane (jako nieporównywalne)

zztt <- left_join(zzt, zz.last.week, by='geo') %>%
  drop_na(lw) %>%
  filter (lw >= latestweek - 6)
## ile zostało krajów (dla ciekawości):
countries.left <- nrow(zztt)

Wykres słupkowy względnej nadmiarowej umieralności (Total excess mortality by country) Trochę go komplikujemy bo słupek PL ma być wyróżniony (w tym celu tworzymy zmienną base a potem ręcznie ustawiamy legendę scale_fill_manual):

p6 <- zztt %>%  mutate( base=ifelse(geo=='PL', "1", "0")) %>% 
  ggplot(aes(x = reorder(country, exmp ), fill=as.factor(base))) +
  geom_bar(aes(y = exmp), stat="identity", alpha=.4 ) +
  geom_text(aes(label=sprintf("%.1f", exmp), y= exmp ), 
            vjust=0.25, hjust=1.25, size=2, color='black' ) +
  geom_text(aes(label=sprintf("%.0f", lw), y= 0 ), 
            vjust=0.25, hjust=-1.25, size=2, color='brown4' ) +
  xlab(label="") +
  ylab(label="") +
  ggtitle("Total Excessive deaths in Europe 2020--2021",
          subtitle ="Sum of (deaths_2020/2021 - average_2015--2019) / average_2015--2019 * 100%") +
  theme_nikw() +
  coord_flip() +
  scale_fill_manual( values = c( "1"="red", "0"="steelblue" ), guide = FALSE ) +
  labs(caption=x.note)

Wykres dynamiki. Ponieważ krajów jest dużo to usuwamy (pierwszy filter) `nieciekawe' (małe lub te które mają jakiś feler, np dane niepełne);

p7 <- zz %>% filter (! geo %in% c('AL', 'AM', 'IE', 'IS', 'UK', 'CY', 'GE', 'LI', 'ME')) %>%
  mutate(date= as.Date(sprintf("%i-%i-1", year, week), "%Y-%U-%u") ) %>%
  ##ggplot(aes(x = date, y =exp, group=geo, color=geo )) +
  ggplot(aes(x = date, y =exp)) +
  facet_wrap( ~country, scales = "fixed", ncol = 4, shrink = F) +
  geom_point(size=.4) +
  geom_smooth(method="loess", size=1, se = F, color="red",  span =.25) +
  geom_hline(yintercept = 50, color="firebrick", alpha=.2, size=1) +
  scale_y_continuous(breaks=seq(-100, 200, by=50)) +
  coord_cartesian(ylim = c(-100, 200)) +
  ggtitle("Excessive deaths in Europe 2020--2021",
          subtitle ="(deaths_2020/2021 - average_2015--2019) / average_2015--2019 * 100%") +
  theme_nikw() +
  xlab(label="%") +
  labs(caption=note)

Teraz wszystko na jednym ale ponieważ krajów jest za dużo (w sensie za mało byłoby kolorów żeby wykres był czytelny), to dzielimy je na dwie grupy: Polska i reszta. Każda grupa innym kolorem:

p8 <- zz %>% filter (! geo %in% c('AL', 'AM', 'IE', 'IS', 'UK', 'CY', 'GE', 'LI', 'ME')) %>%
  mutate( base=ifelse(geo=='PL', "1", "0")) %>% 
  mutate(date= as.Date(sprintf("%i-%i-1", year, week), "%Y-%U-%u") ) %>%
  ggplot(aes(x = date, y =exp,  color=as.factor(base ))) +
  ##ggplot(aes(x = date, y =exp)) +
  geom_point(size=.4) +
  geom_smooth(aes(x = date, y =exp, group=geo), method="loess", size=.3, se = F, span =.25) +
  geom_hline(yintercept = 50, color="firebrick", alpha=.2, size=1) +
  coord_cartesian(ylim = c(-100, 200)) +
  scale_color_manual( values = c( "1"="red", "0"="steelblue" ), guide = FALSE ) +
  theme_nikw() +
  ylab(label="%") +
  ggtitle("Excessive deaths in Europe 2020--2021 and Poland (red)",
          subtitle ="(deaths_2020/2021 - average_2015--2019) / average_2015--2019 * 100%") +
  scale_y_continuous(breaks=seq(-100, 200, by=50)) +
  labs(caption=note)

Wszystko działa (i koliduje jak mówił Bolec), bo działało na PC, ale był problem na raspberry, na którym to ostatecznie ma być uruchamiane z automatu co dwa tygodnie. Bo tam Debian starszy (Buster) i na arma. Zapodanie:

install.package("eurostat")

Skutkowało pobraniem wielu pakietów i błędem przy kompilacji czegoś co się nazywało s2. Kombinowałem na różne sposoby jak to skompilować bo stosownego pakietu w Debianie nie było. Ni-chu-chu.

W przypływie olśnienia obejrzałem zależności tego eurostatu, wśród których żadnego s2 nie było. Był za to sf, które z kolei rzekomo zależał od s2. Żeby było śmieszniej sf dało się zainstalować (mimo że s2 nie było -- stąd rzekomo):

apt install r-cran-sf/oldstable

Problem się rozwiązał...

Skrypt będzie się uruchamiał co dwa tygodnie. Generował rysunki i wysyłał je na twittera

Całość jest na githubie i to w dodatku jako plik Rmd, ale nie na moim koncie tylko na koncie Koła Naukowego SM w PSW w Kwidzynie, którego póki co jestem opiekunem i jedynym członkiem w jednym (por. github.com/knsm-psw/ES-mortality)

url | Fri, 10/12/2021 05:56 | tagi: eurostat, inne, r, statystyka, zgony

Tygodniowe dane nt. zgonów z GUS

GUS się wychylił niespodziewanie z dużą paczką danych nt. zgonów. Dane są tygodniowe, w podziale na płcie, regiony w klasyfikacji NUTS oraz 5 letnie grupy wiekowe.

Dane są udostępnione w formacie XSLX w dość niepraktycznej z punktu widzenia przetwarzania strukturze (kolumny to tygodnie/wiersze to różne kategorie agregacji: płeć, wiek, region), który zamieniłem na CSV o następującej prostej 7 kolumnowej strukturze:

year;sex;week;date;age;geo;value

W miarę oczywiste jest, że year to rok, sex to płeć, week to numer tygodnia, date to data pierwszego dnia tygodnia (poniedziałek), geo to identyfikator obszaru a value liczba zgonów odpowiadająca kolumn 1--6. Ten plik jest podzielony na lata bo w całości zajmuje circa 200Mb. Umieściłem go tutaj.

Skrypt też w R wymodziłem co wizualizuje zgony wg grup wieku oraz województw. Ponieważ kombinacji płeć/wiek/region są setki, moje wykresy dotyczą zgonów ogółem/kobiet/mężczyzn w podziale na grupy wiekowe oraz ogółem w podziale na województwa. Każdy wykres zawiera dwa szeregi: liczbę zgonów w 2020 roku oraz średnią liczbę zgonów z lat 2015--2019. Ponadto jest wykres z jedną krzywą: procent liczony dla stosownych tygodni jako liczba zgonów w 2020 przez średnią 5 letnią z lat 2015--2019. Ten wykres występuje też w wariancie skróconym: tylko 6 ostatnich tygodni, co pozwala dodać do punktów wartości liczbowe (które nie zachodzą na siebie).

library("ggplot2")
library("dplyr")
library("scales")
library("ggthemes")
library("ggpubr")
library("tidyr")

picWd <- 12
spanV <- 0.5
GUS.url <- "https://stat.gov.pl/obszary-tematyczne/ludnosc/ludnosc/zgony-wedlug-tygodni,39,2.html"
NIKW.url <- "(c) NI-KW @ github.com/knsm-psw/GUS_mortality"
NIKW <- sprintf ("%s | %s", GUS, NIKW.url)

z <- read.csv("PL-mortality-2015.csv", sep = ';',  header=T, na.string="NA" )
lastO <- last(z$date)
lastT <- last(z$week)

nuts <- c('PL21', 'PL22', 'PL41', 'PL42', 'PL43', 'PL51', 'PL52', 'PL61', 'PL62',
'PL63', 'PL71', 'PL72', 'PL81', 'PL82', 'PL84', 'PL91', 'PL92')

### Ogółem
z00 <- z %>% filter ( sex == 'O'  & geo == 'PL' ) %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(age,week) %>% summarise( year = 't19',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
zz1 <- z1 %>% group_by(age,week) %>% summarise( year = 't20',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

farbe19 <- '#F8766D'
farbe20 <- '#00BFC4'

p1 <- ggplot(zz1, aes(x=week, y=vv, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~age, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"), values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg grup wiekowych (PL/Ogółem)", subtitle=sprintf("%s | ostatni tydzień: %s", NIKW, lastO) )

ggsave(plot=p1, "zgony_PL_by_age_O.png", width=picWd)

### M ###
z00 <- z %>% filter ( sex == 'M'  & geo == 'PL' ) %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(age,week) %>% summarise( year = 't19',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
zz1 <- z1 %>% group_by(age,week) %>% summarise( year = 't20',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

p2 <- ggplot(zz1, aes(x=week, y=vv, group=year, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~age, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 ##labs(caption=source) +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"),
   values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg grup wiekowych (PL/Mężczyźni)", subtitle=sprintf("%s | ostatni tydzień: %s", NIKW, lastO) )

ggsave(plot=p2, "zgony_PL_by_age_M.png", width=picWd)


### K #########################################
z00 <- z %>% filter ( sex == 'K'  & geo == 'PL' ) %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(age,week) %>% summarise( year = 't19',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
zz1 <- z1 %>% group_by(age,week) %>% summarise( year = 't20',
  vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

p3 <- ggplot(zz1, aes(x=week, y=vv, group=year, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~age, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 ##labs(caption=source) +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"),
   values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg grup wiekowych (PL/Kobiety)", subtitle=sprintf("%s | ostatni tydzień: %s", NIKW, lastO) )

ggsave(plot=p3, "zgony_PL_by_age_K.png", width= picWd)

### ogółem wg województw #####################################
n <- read.csv("nuts.csv", sep = ';',  header=T, na.string="NA" )
## dodaj nazwy
z <- left_join(z, n, by='geo')

## wiek razem
z00 <- z %>% filter ( sex == 'O' & geo %in% nuts & age == 'OGÓŁEM') %>% as.data.frame

z0 <- z00 %>% filter ( year >= 2015  & year < 2020 ) %>% as.data.frame
z1 <- z00 %>% filter ( year == 2020 ) %>% as.data.frame

## średnie w okresie 1 -- (n-1)
zz0 <- z0 %>% group_by(name,week) %>%
 summarise( year = 't19', vv = mean(value, na.rm=TRUE)) %>% as.data.frame
 zz1 <- z1 %>% group_by(name,week) %>%
  summarise( year = 't20', vv = mean(value, na.rm=TRUE)) %>% as.data.frame
### Połącz
zz1 <- bind_rows(zz0, zz1)

lastWeek <- last(zz1$week)
firstWeek <- lastWeek - 6

zz1 <- zz1 %>% filter ( week >= firstWeek  ) %>% as.data.frame
print(zz1)

p4 <- ggplot(zz1, aes(x=week, y=vv, group=year, color=year)) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~name, scales = "free_y") +
 xlab(label="") +
 ylab(label="") +
 ##theme_nikw()+
 ##labs(caption=source) +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok: ", labels = c("średnia 2015--2019", "2020"),
   values = c("t19"=farbe19, "t20"=farbe20 )) +
 ggtitle("Zgony wg województw* (PL/Ogółem)", 
   subtitle=sprintf("*wg klasyfikacji NUTS stąd mazowieckie/stołeczne | %s | ostatni tydzień: %s", NIKW, lastO))

ggsave(plot=p4, "zgony_PL_by_woj_O.png", width=picWd)

## jako %% w średniej w poprzednich 5 lat

zz1 <- zz1 %>% spread(year, vv)

zz1$yy <- zz1$t20 / zz1$t19 * 100

p5 <- ggplot(zz1, aes(x=week, y=yy), color=farbe20) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4, color=farbe20) +
 geom_point(size=.4, alpha=.5) +
 facet_wrap( ~name, scales = "fixed") +
 xlab(label="nr tygodnia") +
 ylab(label="%") +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok 2020: ", labels = c("% 2020/(średnia 2015--2015)"),
   values = c("yy"=farbe20 )  ) +
 ggtitle("Zgony wg województw* (PL/Ogółem)", 
   subtitle=sprintf("*wg klasyfikacji NUTS stąd mazowieckie/stołeczne | %s | ostatni tydzień: %s", NIKW, lastO))

ggsave(plot=p5, "zgony_PL_by_woj_P.png", width=picWd)

zz1 <- zz1 %>% filter ( week >= firstWeek  ) %>% as.data.frame

p6 <- ggplot(zz1, aes(x=week, y=yy), color=farbe20) +
 geom_smooth(method="loess", se=F, span=spanV, size=.4, color=farbe20) +
 geom_point(size=.4, alpha=.5) +
 geom_text(aes(label=sprintf("%.1f", yy)), vjust=-1.25, size=1.5) +
 facet_wrap( ~name, scales = "fixed") +
 xlab(label="nr tygodnia") +
 ylab(label="%") +
 theme(plot.subtitle=element_text(size=9), legend.position="top")+
 scale_color_manual(name="Rok 2020: ", labels = c("% 2020/(średnia 2015--2015)"),
   values = c("yy"=farbe20 )  ) +
   ggtitle(sprintf("Zgony wg województw* (PL/Ogółem) tygodnie: %i--%i (%i tydzień zaczyna się %s)",
     firstWeek, lastWeek, lastWeek, lastO), 
   subtitle=sprintf("*wg klasyfikacji NUTS stąd mazowieckie/stołeczne | %s", NIKW))

ggsave(plot=p6, "zgony_PL_by_woj_P6.png", width=picWd)

url | Mon, 23/11/2020 17:26 | tagi: gus, inne, r, statystyka, zgony

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: eurostat, inne, r, statystyka

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: inne, r, statystyka

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: hhi, inne, mevo, r, statystyka

Hojarskiej ekstremalnie skośny rozkład głosów

Uczono mnie, że współczynnik skośności aczkolwiek może przyjąć wartości większe od 3, to w praktyce taka sytuacja się nie zdarza. Dlatego z rezerwą podszedłem do obliczeń, z których wynikało, że dla pewnego zbioru danych wynosi on 14:

library(moments)
# Wynik 1-ki z listy komitetu Z.Stonogi (D.Hojarska, wybory 2015)
s <- read.csv("hojarska.csv", sep = ';', header=T, na.string="NA");
skewness(s$glosy)
[1] 14.08602

nrow(s)
[1] 657

sum(s$glosy)
[1] 1671

Pierwsza hipoteza jest taka, że formuła liczenia współczynnika może być egzotyczna. Ustalmy zatem jak toto jest liczone:

?skewness
Description:  
This function computes skewness of given data

Wiele to się nie dowiedziałem. Po ściągnięciu źródła można ustalić, że to współczynnik klasyczny czyli iloraz trzeciego momentu centralnego przez odchylenie standardowe do trzeciej potęgi. Zatem jak najbardziej klasyczny a nie egzotyczny. Sprawdźmy co wyliczy Perl dla pewności:

#!/usr/bin/perl

print STDERR "** Użycie: $0 plik numer-kolumny (pierwszy wiersz jest pomijany)\n";
print STDERR "** Domyślnie: wyniki kandydata nr1 na liście komitetu Z.Stonoga (okr.25 / D.Hojarska)\n";

$file = "hojarska.csv"; ## default file
$column = 1; ## first column = 0

if ( $#ARGV >= 0 ) { $file=$ARGV[0]; } 
if ( $#ARGV >= 1 ) { $column=$ARGV[1]; } 

open (S, "$file") || die "Cannot open $file\n";

print "\nDane z pliku $file (kolumna: $column):\n"; 
$hdr = <S> ; # wczytaj i pomin naglowek ($n będzie prawidłowe)
while (<S>) { chomp();
  @tmp = split (/;/, $_); 
  $v = $tmp[$column];
   push(@L, $v) ; 
   $sum += $v; 
   $Counts{"$v"}++;
   $n++;
}

# Wyznaczenie średniej:
$mean = $sum /$n; 

## Wyznaczenie dominanty:
## przy okazji wydrukowanie danych pogrupowanych
print "+--------+---------+--------+--------+--------+\n";
print "| Głosy  | Obwody  |  cumGł | cumGł% |  cumN% |\n";
print "+--------+---------+--------+--------+--------+\n";
$maxc = -1;
for $c (sort {$a <=> $b } keys %Counts ) { 
  $sumCum += $Counts{"$c"} * $c; 
  $nCum += $Counts{"$c"}; 
  if ($maxc_pos < $Counts{$c} ) { 
    $maxc_pos = $Counts{$c} ; $maxc_val = $c
  }
  printf "| %6d |  %6d | %6d | %6.2f | %6.2f |\n", 
    $c,  $Counts{$c}, $sumCum, $sumCum/$sum *100, $nCum/$n *100;
}
print "+--------+---------+--------+--------+--------+\n\n";

$mode = $maxc_val; # dominanta

## Wyznaczenie mediany:
$half = int(($#L +1)/2 );
@L = sort ({ $a <=> $b } @L) ; #numerycznie
##print "$half of @L\n";
$median = $L[$half];

printf "Średnie: x̄ = %.3f (N=%d) Me =%.3f D = %.3f\n", $mean, $n, $median, $mode;

## Odchylenia od średniej:
for my $l (@L) {##
  $sd1 = ($l - $mean) ; 
  $sd2 = $sd1 * $sd1; # ^2 
  $sd3 = $sd2 * $sd1; # ^3

  $sum_sd1 += $sd1; $sum_sd2 += $sd2;
  $sum_sd3 += $sd3;
}

# odchylenie std./3-ci moment:
$sd = sqrt(($sum_sd2 / $n)); 
$u3 = $sum_sd3/$n;

printf "Sumy (x-x̄): %.2f(¹) %.2f(²) %.2f(³)\n", 
   $sum_sd1, $sum_sd2, $sum_sd3;
printf "Rozproszenie: σ = %.3f µ³ = %.3f\n", $sd, $u3;

printf "Skośność: (x̄-D)/σ = %.2f", ($mean -$mode)/$sd;
printf " 3(x̄-Me)/σ) = %.2f", 3 * ($mean - $median)/$sd;
printf " (µ³/σ³) = %.2f\n", $u3/($sd*$sd*$sd);
##//

Uruchomienie powyższego skryptu daje w wyniku to samo. Moja pewność, że wszystko jest OK jest teraz (prawie że) stuprocentowa.

Wracając do R:

library(ggplot2);
p <-  ggplot(data = s, aes(x = glosy)) + geom_histogram(binwidth = 4)
p
breaks <- ggplot_build(p)$data
breaks

[[1]]
     y count   x xmin xmax      density      ncount    ndensity PANEL group
1  482   482   0   -2    2 0.1834094368 1.000000000 2628.000000     1    -1
2  140   140   4    2    6 0.0532724505 0.290456432  763.319502     1    -1
...

Pierwszy przedział jest definiowany jako (-2;2], ale z uwagi na wartość danych de facto liczy głosy w komisjach, w których Hojarska dostała 0--2 głosów (drugi to (2;6], itd) i ni cholery nie da się tego zmienić na coś bardziej zbliżonego do prawdy. Bo tak jak jest to faktycznie pierwszy przedział jest dwa razy węższy niż każdy pozostały. W szczególności prawdziwa średnia, w tym przedziale wynosi 1,259 głosa, zaś liczona jako środek przedziału przez liczebność wynosi 1.0 (482 *1/482). Wg formuły ggplota środek przedziału to zero zatem średnia też wyjdzie 0, tj. wartość jest wyznaczona z dużym/nieokreślonym błędem.

Dzieje się tak ponieważ: the histogram is centered at 0, and the first bars xlimits are at 0.5*binwidth and -0.5*binwidth. Dopiero gdy dane są dodatnie, to ggplot zaczyna od zera. No ale nasz zbiór zawiera zera i klops. Zmiana tego jest nietrywialna.

Zamiast ggplota moża użyć po prostu polecenia hist:

h <- hist(s$glosy, breaks=seq(0,max(s$glosy), by=4) )
h

$breaks
 [1]   0   4   8  12  16  20  24  28  32  36  40  44  48  52  56  60  64  68  72
[20]  76  80  84  88  92  96 100 104 108 112 116 120 124 128 132 136 140 144 148
[39] 152

$counts
 [1] 592  42   6   4   2   2   2   1   2   1   2   0   0   0   0   0   0   0   0
[20]   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1

$density
 [1] 0.2252663623 0.0159817352 0.0022831050 0.0015220700 0.0007610350
 [6] 0.0007610350 0.0007610350 0.0003805175 0.0007610350 0.0003805175
[11] 0.0007610350 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[16] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[21] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[26] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[31] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[36] 0.0000000000 0.0000000000 0.0003805175

$mids
 [1]   2   6  10  14  18  22  26  30  34  38  42  46  50  54  58  62  66  70  74
[20]  78  82  86  90  94  98 102 106 110 114 118 122 126 130 134 138 142 146 150

No teraz liczy (prawie) jak trzeba pierwszy przedział [0;4], drugi (4;8] itd. Prawie bo na upartego pierwszy przedział jest szeroki na 5 wartości a kolejne na 4. Eh te detale i te upierdliwe zero:-)

Przy okazji się zadumałem, a cóż to jest to density. W pierwszym podejściu, że to są częstości, ale nie, bo ewidentnie nie sumują się do 1:

  ?hist
  density: values f^(x[i]), as estimated density values. If
          `all(diff(breaks) == 1)', they are the relative frequencies
          `counts/n' and in general satisfy sum[i; f^(x[i])
          (b[i+1]-b[i])] = 1, where b[i] = `breaks[i]

No dość kryptycznie, ale można się domyśleć, że nie $\sum_i d_i =1$, ale $\sum_i d_i * w_i$, gdzie $w_i$, to szerokość przedziału $i$ (a $d_i$ oznacza gęstość w przedziale $i$ oczywiście). Pole pod krzywą będzie zatem równe 1, jak tego należy się spodziewać po gęstości.

Wprawdzie już to padło mimochodem wyżej, ale analizowany zbiór danych to liczba oddanych głosów na numer jeden na liście komitetu wyborczego Zbigniewa Stonogi w okręgu wyborczym 25 w wyborach do Sejmu w 2015 roku. Tym numerem jeden była Danuta Hojarska, która faktycznie dostała 152 głosy w jednej komisji (we wsi, w której mieszka) oraz 0--4 głosów w 90% pozostałych komisjach (w tym 0 głosów w 191/657.0 = 29% komisjach).

Analiza przestrzenna

Ze starych zapasów odgrzebałem plik CSV ze współrzędnymi obwodowych komisji wyborczych i połączyłem go z wynikami uzyskanymi przez Hojarską z zamiarem wyświetlenia tego na mapie za pomocą Google Fusion Tables (GFT):

obwod;coordexact;glosy;icon
101884;54.2701205 18.6317438;3;measle_white
101885;54.260714 18.6282809;0;measle_white
101886;54.262187 18.632082;2;measle_white
101887;54.257501 18.63147;1;measle_white
101888;54.25786 18.6306574;7;measle_grey
...

Kolumna icon zawiera nazwę ikony, która jest definiowana wg schematu: 0--3 głosy biala (measle_white); 4--9 głosów szara (measle_grey) 9-19 głosów żółta (small_yellow) 20--99 głosów czerwona (small_red) 100 i więcej głosów purpurowa (small_purple). Zestawienie dostępnych w GFT ikon można znaleźć tutaj. Konfigurowanie GFT aby korzystała z kolumny z nazwami ikon sprowadza się do wyklikania stosownej pozycji (Configure map →Feature map →Change feature styles →Use icon specified in a column.)

Przy okazji odkurzyłem zasób skryptów/zasobów używanych do geokodowania i/lub przygotowywania danych pod potrzeby GFT, w szczególności: joincvs.pl (łączy dwa pliki csv w oparciu o wartości n-tej kolumny w pierwszym pliku oraz m-tej kolumny w drugim); geocodeCoder0.pl (uruchamia geokodera Google). Oba skrypty i jeszcze parę rzeczy można znaleźć: tutaj.

url | Sun, 25/02/2018 17:38 | tagi: gft, statystyka

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: r, statystyka, twitter

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: inne, r, statistics, statystyka

Numerologia stosowana (lies, damned lies and statistics)

,,Po dziennik Fakt (Ringier Axel Springer Polska) w okresie od maja do października 2015 roku sięgało 10,33%. Drugie miejsce zajęła Gazeta Wyborcza (Agora), której poziom czytelnictwa wzrósł w ciągu roku o 0,39 pkt. proc. do 8,58%.'' (za Gazeta Wyborcza dogania...).

Na samym końcu tekstu: ,,Podane liczby to wskaźniki CCS -- Czytelnictwo Cyklu Sezonowego''.

Ki ch#j ten CCS? Bo liczby wydają się mooocno napompowane. Okazuje się, że CCS oznacza jaki procent w populacji stanowią osoby, które zetknęły się z tym tytułem choćby raz w okresie cyklu sezonowego... (za Wskaźniki PBC). Jeszcze żeby było wiadomo co oznacza ,,zetknęły się''? Zawinęły karpia w sklepie albo pocięły jako ersatz Toilettenpapier?

Skądinąd wiadomo (z innej strony na tym samym serwerze), że sprzedaż ww. jest licha (por. Spadła sprzedaż kioskowa wszystkich dziennikóww). Proste przemożenie 30136052 (Polacy w wieku 15--75 lat) x 10,33% i 8,58% odpowiednio daje w wyniku 3,1 i 2,6 mln czytelników przy 300 i 180 tys sprzedanych egz. co oznacza, że 1 egz. Faktu czyta 10 osób, a GW nawet 15... Kto chce niech wierzy.

Ja bym zamiast pytać ,,czy się zetknął'' pytał się o treść 1 dowolnego tekstu z gazety ,,z którą się miało zetknięcie''... Ciekawe ile by wtedy wyszło...

Dane z artykułów w formacie CSV są poniżej:

tytul;CCS2015;CCS2014;platne2015;druk2015;platne2014;druk2014;cz2015;cz2014;razem2015u;razem2014u;druk2015u;druk2014u
Fakt;10.33;11.74;308616;305178;326748;323002;3113054;3537973;1008.7;1082.8;1020.1;1095.3
Gazeta Wyborcza;8.58;8.19;177447;106794;191127;113292;2585673;2468143;1457.2;1291.4;2421.2;2178.6
SuperExpress;3.58;4.55;145647;143954;153387;151735;1078871;1371190;740.7;893.9;749.5;903.7
Przegląd Sportowy;3.57;2.90;31484;30518;35546;34500;1075857;873946;3417.2;2458.6;3525.3;2533.2
Dziennik Gazeta Prawna;1.68;1.59;54089;9887;58027;11412;506286;479163;936.0;825.8;5120.7;4198.8
Rzeczpospolita;1.62;1.94;56560;7897;58056;9352;488204;584639;863.2;1007.0;6182.1;6251.5
PulsBiznesu;0.31;0.20;13314;1655;13996;1743;93422;60272;701.7;430.6;5644.8;3458.0
Parkiet Gazeta Giełdy;0.09;0.08;5332;1319;5015;1532;27122;24109;508.7;480.7;2056.3;1573.7

Gdzie: CCS2015 -- wskaźniki CCS dla 2015 r.; platne2015 -- rozpowszechnianie płatne razem w 2015; druk2015 -- sprzedaż egz. wydań drukowanych; cz2015 = CCS2015 * 30136052; razem2015u = platne2015 / druk2015 * 100%; druk2015u = druk2015 / druk2015 * 100%.
Dla roku 2014 analogicznie.

Ponadto:

Sprzedaż ogółem to suma sprzedaży egzemplarzowej wydań drukowanych, sprzedaży egzemplarzowej e-wydań oraz wszystkich form prenumeraty wydań drukowanych i prenumeraty e-wydań.

Rozpowszechnianie płatne razem to suma sprzedaży egzemplarzowej wydań drukowanych, sprzedaży egzemplarzowej e-wydań, prenumeraty wydań drukowanych i prenumeraty e-wydań oraz innych płatnych formy rozpowszechniania wydań drukowanych i innej płatnej dystrybucji e-wydań.

url | Mon, 28/12/2015 20:03 | tagi: inne, statistics, statystyka

Afera madrycka: taka tam analiza wyjazdów posłów 7 kadencji

UWAGA: Ten tekst nie jest o polityce ale o [elementarnej] statystyce.

Media informowały, że posłowie PiS Adam Hofman, Mariusz A. Kamiński i Adam Rogacki wzięli na podróż do Madrytu na posiedzenie komisji Zgromadzenia Parlamentarnego Rady Europy po kilkanaście tysięcy złotych zaliczki, zgłaszając wyjazd samochodem; w rzeczywistości polecieli tanimi liniami lotniczymi. Ponieważ kontrola wydatków posłów jest iluzoryczna różnica pomiędzy kosztem podróży samochodem a samolotem [za dużo mniejsze pieniądze] miała stanowić dodatkowy przychód wyżej wymienionych. Według prokuratury, która wszczęła śledztwo, zachodzi podejrzenie popełnienia oszustwa.

Łapiąc wiatr w żagle [sprawa się upubliczniła tuż przed ostatnimi wyborami samorządowymi] koalicja rządząca w osobie Marszałka Sejmu RP Sikorskiego zarządziła audyt, którego efektem było udostępnienie m.in. dokumentu pn. Wyjazdy zagraniczne posłów VII kadencja (kopia jest tutaj).

Jak przystało na kraj, w którym od lat działa Ministerstwo cyfryzacji zestawienie jest w formacie PDF, zatem pierwszym ruchem była zamiana na coś przetwarzalnego. Wpisanie w google PDF+Excel+conversion skutkuje ogromną listą potencjalnych konwerterów. Bagatelizując skalę problemu spróbowałem dokonać konwersji narzędziami dostępnymi on-line, ale z marnym rezultatem (za duży dokument przykładowo; serwis za free jest tylko dla PDFów mniejszych niż 50 stron). W przypadku Convert PDF to EXCEL online & free coś tam skonwertował, nawet wyglądało toto na pierwszy rzut oka OK ale na drugi już nie: dokument niekompletny oraz nieprawidłowo zamienione niektóre liczby (przykładowo zamiast 837,50 zł w arkuszu jest 83750 -- 100 razy więcej!).

Ostatecznie skończyło się na ściągnięciu 30 dniowej wersji Adobe Acrobata Pro XI, który faktycznie sprawdził się w roli konwertera PDF→XLSX. Do konwersji wykorzystałem służbowego laptopa Elki wyposażonego w legalny Office 2010, na którym zainstalowałem ww. AA Pro XI. OOffice niby czyta XLSX, ale z koszmarnymi błędami, więc żeby dalej móc obrabiać arkusz w Linuksie wczytałem wynikowy XLSX do Excela 2010 po czym zapisałem go w (starszym) formacie XLS. Ten plik wyświetlił się w OO Calcu bez problemu.

Arkusz jest tak sformatowany, że 4 pierwsze komórki oraz są często wielowierszowe i scalone, zawierają bowiem liczbę porządkową, datę, miejsce i cel wyjazdu delegacji posłów. Po zamianie na plik CSV zawartość komórek scalonych pojawi się w pierwszym wierszu, a pozostałe będą puste. Prostym skryptem Perlowym mogę wypełnić puste komórki wg. algorytmu: jeżeli cztery pierwsze pola są puste, to skopiuj wartości ostatnich niepustych:

if ($tmp[0] eq '' && $tmp[1] eq '' && $tmp[2] eq '' && $tmp[3] eq '' ) { ... }

Pierwszy problem: wielowierszowe komórki z kolumn 1--4 nie zawsze są scalone. Czasem tekst jest podzielony na wiersze co psuje konwersję. Ręcznie scalam niescalone komórki (trochę to trwa). Przed scaleniem usuwam z kolumn 1--4 końce wiersza.

Drugi problem: część liczb nie jest liczbami z uwagi na użycie separatora tysięcy, który się zamienił w PDFie na odstęp (spację). Zatem zaznaczam kolumny zawierające różne pozycje kosztów po czym:

Edytuj→Znajdź i zamień: usuwam odstępy, tj. zamieniam spację na pusty napis
Format→Komórki: wybieram numer z dwoma miejscami po przecinku.

Po uporządkowaniu arkusza, zapisuję go w formacie CSV. Następnie prostym skryptem Perlowym zamieniam na taki plik CSV, w którym puste komórki są wypełniane zawartością z poprzednich wierszy. Kolumna Państwo - miasto jest kopiowana. Kopia jest zmieniana na jednoznaczne: Państwo, miasto (pierwszy-kraj, przecinek, pierwsze miasto z listy celów podróży -- żeby geokoderowi było łatwiej.)

Innym skryptem Perlowym dodaję do pliku CSV 3 kolumny, które zawierają:

współrzędne celu podróży (w tym celu zamieniam adres Państwo, miasto na współrzędne geograficzne korzystając z geokodera Google);
odległość w kilometrach pomiędzy punktem o współrzędnych 21.028075/52.225208 (W-wa, Wiejska 1) a celem podróży (obliczoną przy wykorzystaniu pakietu GIS::Distance);
linię zdefiniowana w formacie KML o końcach 21.028075/52.225208--współrzędne-celu-podróży (do ewentualnego wykorzystania z Google Fusion Tables).

#!/usr/bin/perl
#
use Storable;
use Google::GeoCoder::Smart;
use GIS::Distance;

$geo = Google::GeoCoder::Smart->new();

my $gis = GIS::Distance->new();

my $GeoCodeCacheName = 'geocode.cache';
my $NewCoordinatesFetched=0; # global flag
my $SLEEP_TIME = 2 ;
my $coords_okr = "21.028075,52.225208"; # Warszawa = środek świata

my %GeoCodeCache = %{ retrieve("$GeoCodeCacheName") } if ( -f "$GeoCodeCacheName" ) ;
my ($wwa_lng, $wwa_lat) = split (",", $coords_okr);
my $linesNo = 0 ;
my $GCtotaluse = 1; # laczna liczba wywolan geocodera

while (<>) {
  $linesNo++;
  chomp();  $_ =~ s/[ \t]+;[ \t]+/;/g; ## usuń ew. niepotrzebne spacje

  @line = split ";", $_;  print STDERR "**$linesNo = $line[3] ... ";

  # geokodowanie (uwaga na limit) 
  # Poprawki dla miejsc, których nie zna Google:
  $line[3] =~ s/Erewań/Erywań/; ## 
  $line[3] =~ s/Sowayma/Madaba/; ## najbliższe miasto
  $line[3] =~ s/Bołszowce/Iwano-Frankiwsk/; ## najbliższe miasto

  my $coords = addr2coords( $line[3] );

  ($tmp_lat, $tmp_lng, $gcuse) = split " ", $coords;
  if ($gcuse > 0) {$GCtotaluse++ ; }

  $distance = $gis->distance($tmp_lat,$tmp_lng => $wwa_lat,$wwa_lng );
  $distance_t = sprintf ("%.1f", $distance);

  my $kml_line = "<LineString><coordinates>$tmp_lng,$tmp_lat $coords_okr</coordinates></LineString>";
  print "$_;\"$coords\";$distance_t;\"$kml_line\"\n";
  print STDERR "\n";

  if ($GCtotaluse % 100 == 0 ) {# store every 100 geocoder calls
    store(\%GeoCodeCache, "$GeoCodeCacheName");
    print STDERR "\n... Cache stored. ***\n";    
  }
}

##
store(\%GeoCodeCache, "$GeoCodeCacheName");

## ## ## ####
sub addr2coords {
 my $a = shift ;
 my $r = shift || 'n';
 my ($lat, $lng) ;
 my $GCuse = 0;

 ##consult cache first
 if (exists $GeoCodeCache{"$a"} ) {
   print STDERR "Coordinates catched ... $a ";
   ($lat,$lng) = split (" ", $GeoCodeCache{"$a"} );
 }
 else {
   print STDERR "Geocoding ... $a ";
   my ($resultnum, $error, @results, $returncontent) = $geo->geocode("address" => "$a");
   $GCuse = 1;
   sleep $SLEEP_TIME; ## make short pause

   $resultnum--; 
   $resultNo=$resultnum ;

   if (resultNo > 0) { print STDERR "** Location $a occured more than once! **" }
   if ($error eq 'OK') {
     $NewCoordinatesFetched=1;
     for $num(0 .. $resultnum) {
       $lat = $results[$num]{geometry}{location}{lat};
       $lng = $results[$num]{geometry}{location}{lng};
       ##print "*** LAT/LNG:$lat $lng ERROR: $error RES: $resultNo ***\n";
     }

     $GeoCodeCache{"$a"} = "$lat $lng"; ## store in cache

   } else { print STDERR "** Location $a not found! due to $error **"  }
 }


 if ($r eq 'r' ) { return "$lng,$lat,$GCuse"; } # w formacie KML
 else { return "$lat $lng $GCuse"; }
}

Gotowy plik CSV zawierający zestawienie podróży jest dostępny tutaj.

Na podstawie zestawienia i z użyciem pakietu ggplot2 generują się takie oto śliczne wykresy.

Wszystkie podróże z zestawienie (N=1874; odpowiednio: koszt łączny, koszt transportu, długość w tys km):

Tylko podróże dla których koszt transportu był niezerowy (N=1423; odpowiednio: koszt łączny, koszt transportu, długość w tys km):

Poniższy skrypt R sumuje i drukuje wszystkie podróże każdego posła:

require(plyr)

d <- read.csv("W7RR_podroze_by_podroz1.csv", sep = ';', dec = ",",  header=T, na.string="NA");

# Dodaj kolumnę której wartości to konkatenacja: "Poseł|Klub"
d[,"PosKlub"] <- do.call(paste, c(d[c("Posel", "Klub")], sep = "|"));

# Usuń wszystko za wyjątkiem tego co potrzeba:
d <- d[ c("PosKlub", "Klacznie", "Ktransp", "Dist") ];

# Sumowanie po PosKlub 
PSums <- as.data.frame ( ddply(d, .(PosKlub), numcolwise(sum)) );

# Z powrotem rozdziel kolumnę "Poseł|Klub" na dwie
PSums <- as.data.frame ( within(PSums, PosKlub <-data.frame( do.call('rbind', 
   strsplit(as.character(PosKlub), '|', fixed=TRUE))))  )

# Drukuj 
PSums;

Z pliku .Rout kopiuję zestawienie łącznych wydatków posłów oraz łącznej pokonanej przez nich odległości:

       PosKlub.X1 PosKlub.X2 KlacznieT  KtranspT    DistT
1 Adam Abramowicz        PiS   4.02599   2.64595   1.3153
2     Adam Hofman        PiS 119.55271  59.53315  26.1716
3   Adam Kępiński        SLD  10.15754   7.93882   3.8069
4   Adam Kępiński         TR  12.63098   8.02327   2.2107
...

Uwaga: kilkanaście nazwisk się powtarza ponieważ posłowie zmienili przynależność klubową w czasie trwania kadencji [Aby uwzględnić takich posłów sumowanie odbywało się po wartościach zmiennej zawierającej połączone napisy Poseł|Klub.]

Na podstawie takiego z kolei zestawienia i znowu z użyciem ggplot2 generują inne śliczne wykresy.

Uwaga: sumowane tylko podróże, dla których koszt transportu był niezerowy (N=1423; odpowiednio: koszt łączny, koszt transportu, długość w tys km):

Link do tabeli zawierającej zestawienie podróży w formacie Google Fusion Tables jest tutaj.

Dane + skrypty dostępne są także w: github.com/hrpunio/Data.

url | Tue, 09/12/2014 19:09 | tagi: gft, google, inne, perl, politics, polityka, statistics, statystyka

Data sets for statistical education

IMHO there is acute lack of interesting data sets for statistical education purposes, particularly cross-sectional in nature and not related to biology/medicine/astronomy. The domains I am looking for are management, demography and economics but sport is OK too as many students are interested in sport.

I agree there is a lot of data published (mainly) by various government agencies but its almost 100% time series data. There are many test datasets too (like UC Irvine Machine Learning Repository), but they are usually boring, often categorical and short:-)

To remedy a little the lack of good data I compiled myself a few databases and I intend to add more in the future. At present my repository consists of 4 datasets concerning: -- 1100 or so German WW2 U-boots; -- 2200 or so Titanic passangers/crew members -- 1200 or so Rugby Union players who took part in RWC 2011 and RWC 2007 tournaments.

The repository is available at my github account.

url | Tue, 08/11/2011 19:52 | tagi: inne, r, statistics, statystyka

Średnia ruchoma aka moving average

Poniżej nie do końca sprawdzona procedura wyznaczania średniej ruchomej zwykłej, zaimplementowana w Perlu:

#!/usr/bin/perl
my @qq= (1, 2, 6, 7, 10, 11, 11, 12, 13); # 9 elementów

print "Oryginalne:: @qq\n";

$ref_av = moving_avg(\@qq, 3);
print "Uśrednione:: @{$ref_av}\n";

$ref_av = moving_avg(\@qq, 7);
print "Uśrednione:: @{$ref_av}\n";

# ## ## ## ## ## ## ##
sub moving_avg {
  my $wlista = shift; # wskaźnik do listy, która będzie uśredniona
  my $lokr = shift; # liczba okresów

  my $i; my $total ;

  my @lista_av = @{$wlista};

  for ($i=0; $i < $lokr; $i++) {  $total += ${$wlista}[$i] ; }

  ## trzeba sprawdzić czy jest OK; nie jestem pewien
  for ($i = $lokr; $i <= $#{$wlista} + 1; $i++) {
    $lista_av[$i] = $total / $lokr ;
    $total += ${$wlista}[$i] - ${$wlista}[$i - $lokr];
  }

  return \@lista_av; ## zwróć wskaźnik do listy z uśredn. wartościami 
}

url | Thu, 10/02/2011 19:54 | tagi: inne, perl, statystyka

Introduction to Data Technologies

Przypadkowo znalazłem całkiem grubą książkę nt. Introduction to Data Technologies. Na super rewelację toto nie wygląda, ale w wolnej chwili obejrzę dokładniej...

url | Thu, 25/11/2010 20:13 | tagi: r, statystyka

Rysowanie diagramów ścieżkowych za pomocą MetaPosta

Diagram ścieżkowy (path diagram)-- popularna notacja służąca do graficznego definiowania równań strukturalnych -- to rodzaj grafu, w którym kółka i prostokąty oznaczają węzły a strzałki krawędzie. Do tego zarówno kółka/prostokąty jaki i strzałki mogą być oznaczone symbolami. No i tu pech, bo tradycyjnie używa się greckich liter, do tego z subskryptami.

Dia to fajny program, ale przy bardziej zaawansowanych diagramach wychodzą jego ograniczenia, jak przykładowo niemożność umieszczania formuł matematycznych. Wstawianie symboli greckich próbowałem obejść modyfikując LOCALE i klawiaturę (da się) ale subskrypty/superskrypty to już się okazały nie do przejścia. Kombinowałem zatem jakby tu łatwo i przyjemnie zrobić diagramy ścieżkowe w czymś innym. W szczególności książka Skrondala i Rabe-Hesketh (Generalized latent variable modeling: multilevel, longitudinal, and structural equation models, ISBN: 9781584880004) -- ewidentnie składania LaTeXem -- zawiera bardzo schludne diagramy ścieżkowe. Podobno też w R coś tam się samo składa... Nic wszakże konkretnego nie ustaliłem i stanęło jak zwykle na MetaPoście.

Tutaj umieszczone bardzo proste makra pozwalają na rysowania diagramów ścieżkowych. Przykładowo narysowanie jednoczynnikowego modelu pomiaru, zawierającego trzy miary refleksyjne będzie wyglądało jakoś tak:

z1  = (40mm,20mm); %% srodek koła
%% Rysowanie czynnika eta_1 z błędem delta_1
ErCircle(z1, "$\eta_1$", "$\delta_1$", "r");

%% miary (refleksyjne):
z15 = (10mm,0mm);  
%% miara x_1 z błędem \epsilon_1
ErRectangle(z15, "$x_1$", "$\epsilon_1$", "l"); 
CircleToRectangle(z15,z1, "$\lambda_1$");

z10 = (10mm,10mm); 
%% miara x_2 z błędem \epsilon_2
ErRectangle(z10, "$x_2$", "$\epsilon_2$", "l"); 
CircleToRectangle(z10,z1, "$\lambda_2$");

z11 = (10mm,20mm); 
%% miara x_3 z błędem \epsilon_3
ErRectangle(z11, "$x_3$", "$\epsilon_3$", "l"); 
CircleToRectangle(z11,z1, "$\lambda_3$");

Bardziej zaawansowany model jest obok na rysunku. Może komuś się też się przyda...

url | Sun, 19/07/2009 00:34 | tagi: inne, sem, statystyka