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.
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)
Darrell Huff. How to lie with statistics (142 strony/a5) https://en.wikipedia.org/wiki/How_to_Lie_with_Statistics
Ż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)
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)
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)
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)")
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...
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).
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.
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)
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.)
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
,,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ń.
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:
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.
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.
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 }
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...
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...