Przygotowuję trzecie wydanie ,,Przewodnika po pakiecie R”. Zmiany w stosunku do wydania drugiego są spore, kilka rzeczy zostało usuniętych, kilka dodanych.
Jedna z dodanych rzeczy to bardziej rozbudowany rozdział o funkcjach z rodziny apply i plyr. Poniżej przeklejam część nowego wydania ,,Przewodnika …” poświęconą pakietowi plyr. Przepraszam też za ewentualne usterki związane z konwersją LaTeXa do HTMLa.
Fragment książki ,,Przewodnik po pakiecie R”, Przemysław Biecek, wydawnictwo GiS, wydanie 3.
W podejściu z użyciem pary plyr+reshape2, pakiet reshape2 służy to zmiany struktury zbioru danych, a pakiet plyr służy do automatyzacji obliczeń, na odpowiednio przygotowanej strukturze danych. Poniżej przedstawiam krótkie wprowadzenie. Osobom zainteresowanym poznaniem dalszych szczegółów polecam prace http://www.jstatsoft.org/v21/i12/paper i http://www.jstatsoft.org/v40/i01/paper.
Zacznijmy od przedstawienia pakietu reshape2. W tym pakiecie godne uwagi są trzy funkcje: reshape2::melt(), reshape2::acast() i reshape2::dcast(). Funkcja melt() ,,roztapia” dane z postaci tabelarycznej (również dla tabel wielowymiarowych) do tzw. postaci wąskiej (nazywanej też postacią rzadką lub postacią wiersz-kolumna-wartość (RCV od ang. row column value). Funkcja acast() pozwala na przekształcenie danych z postaci wąskiej na postać tabeli krzyżowej (przestawnej), gdzie wartościami pól tej tabeli są wartości zadanej funkcji wyznaczonej na podzbiorze danych.
Przedstawmy obie funkcje na przykładzie. Wykorzystamy pakiet SmarterPoland do pobrania z internetowych baz Eurostatu danych o częstości wypadków drogowych w różnych latach, różnych krajach. Takie dane znajdują się w tabeli tsdtr420 w bazie Eurostatu.
Poniżej wykorzystujemy funkcję SmarterPoland::getEurostatRCV() do pobrania danych z tabeli tsdtr420.
library(SmarterPoland) wypadkiWaska <- getEurostatRCV("tsdtr420") ## trying URL 'http://epp.eurostat.ec.europa.eu/NavTree_prod/everybody/BulkDownloadListing?sort=1' file=data%2Ftsdtr420.tsv.gz' ## Content type 'application/x-gzip' length 2355 bytes ## opened URL ## ================================================== ## downloaded 2355 bytes |
Kilka pierwszych wierszy z pobranych danych. W kolejnych kolumnach są: nazwa opisywanego współczynnika, identyfikator kraju, roku oraz liczba ofiar wypadków.
head(wypadkiWaska, 4) ## victim geo time value ## 1 KIL AT 1991 1551 ## 2 KIL BE 1991 1873 ## 3 KIL BG 1991 1114 ## 4 KIL CY 1991 103 |
Podsumowanie danych z częstościami wystąpień zmiennych czynnikowych.
summary(wypadkiWaska) ## victim geo time value ## KIL :532 AT : 38 1991 : 56 Min. : 4.0 ## KIL_MIO_POP:532 BE : 38 1992 : 56 1st Qu.: 117.0 ## BG : 38 1993 : 56 Median : 204.5 ## CY : 38 1994 : 56 Mean : 2070.9 ## CZ : 38 1995 : 56 3rd Qu.: 1066.2 ## DE : 38 1996 : 56 Max. :75426.0 ## (Other):836 (Other):728 NA's :28 |
Zmienna wypadkiWaska jest już w postaci roztopionej. Pierwsze trzy kolumny opisują trzy różne wymiary. Tutaj są to: rodzaj statystyki (liczba ofiar i liczba ofiar przeliczona na milion mieszkańców), identyfikator kraju, rok, którego dotyczy statystyka. Czwarta kolumna określa wartość zadanej statystyki dla wybranego kraju w wybranym roku. W tym przypadku postać wąska opisuje trzy wymiary zmiennych, ale można rozważać zarówno dane o większej jak i o mniejszej liczbie wymiarów.
Wykorzystajmy teraz funkcję dcast() aby przekształcić te dane do postaci tabelarycznej z krajami w wierszach i latami w kolumnach. Przyjmijmy, że interesuje nas wyłącznie statystyka KIL_MIO_POP (,,ofiary na milion mieszkańców”).
wypadkiWaska <- subset(wypadkiWaska, victim=="KIL_MIO_POP") wypadkiSzeroka <- dcast(wypadkiWaska, geo ~ time, mean, na.rm=TRUE) |
Pierwszym argumentem jest ramka danych w postaci roztopionej. Drugim jest formuła opisująca, które wymiary chcemy by były w wierszach (lewa strona formuły) a które w kolumnach (prawa strona formuły). Argument subset określa, które wiersze z wejściowego zbioru danych brać pod uwagę. Trzeci argument wskazuje funkcję, która będzie zastosowana do wszystkich wierszy ramki wypadkiWaska wskazujących na ten sam rok i kraj. W rozważanym przypadku rok i kraj wyznaczają dokładnie jeden wiersz, więc użycie funkcji mean() może wydawać się nadmiarowe, za każdym razem jest to bowiem liczenie średniej z jednej wartości.
dim(wypadkiSzeroka) ## [1] 28 20 wypadkiSzeroka[wypadkiSzeroka$geo %in% c("UK", "SK", "FR", "PL", "ES", "PT", "LV"),1:10] ## geo 1991 1992 1993 1994 1995 1996 1997 1998 1999 ## 10 ES 227 200 163 143 146 139 142 150 144 ## 13 FR 184 173 172 157 154 147 145 153 145 ## 19 LV 375 298 280 305 264 241 232 280 272 ## 22 PL 207 181 165 175 179 165 189 183 174 ## 23 PT 323 310 271 251 271 272 250 210 200 ## 27 SK 116 128 110 119 123 115 146 152 120 ## 28 UK 83 76 69 66 65 64 64 61 61 |
Dla ćwiczeń roztopimy teraz zbiór danych wypadkiSzeroka, sprowadzając go z powrotem do postaci wąskiej z użyciem funkcji melt().
wypadkiRoztopiona <- melt(wypadkiSzeroka, id = "geo", measure=paste(1991:2008)) head(wypadkiRoztopiona) ## geo value time ## X1991 AT 201 1991 ## X1991.1 BE 188 1991 ## X1991.2 BG 129 1991 ## X1991.3 CY 175 1991 ## X1991.4 CZ 129 1991 ## X1991.5 DE 142 1991 |
Teraz omówmy wybrane, przydatne funkcje z pakietu plyr. Tabela [plyrTabela] przedstawia nazwy wszystkich interesujących funkcji. Jednak nie będziemy ich wszystkich omawiać, ponieważ ich sposób działania jest bardzo podobny i wystarczy omówić kilka wybranych przykładów by wiedzieć jak korzystać z pozostałych.
| wejście \ wyjście | macierz | ramka danych | lista | nic |
| macierz | aaply() | adply() | alply() | a\_ply() |
| ramka danych | daply() | ddply() | dlply() | d\_ply() |
| lista | laply() | ldply() | llply() | l\_ply() |
| n powtórzeń | raply() | rdply() | rlply() | r\_ply() |
| zbiór argumentów | maply() | mdply() | mlply() | m\_ply() |
Tabela: Wybrane funkcje z~pakietu \texttt{plyr} do przetwarzania potokowego. W~kolumnach zaznaczono jakiego typu jest wynik funkcji, czy jest to macierz, ramka danych, lista czy też funkcja nie zwraca wartości. W~wierszach zaznaczono co jest argumentem wejściowym, czy jest to macierz danych, ramka danych, lista zbiorów~danych.
Schemat działania funkcji z tego pakietu jest określany ,,dziel/przekształć/łącz” (ang. split/apply/combine). Przetwarzanie potokowe oznacza, że ta sama funkcja jest stosowana do różnych podzbiorów wejściowego zbioru danych. W pierwszym kroku ten zbiór danych jest dzielony na podzbioru najczęściej z uwagi na wybrane zmienne grupujące (etap ang. split), następnie do każdego wydzielonego podzbioru stosowana jest zadana funkcja (etap ang. apply) w ostatnim kroku otrzymane wyniki są łączone (etap ang. combine).
Konwencja nazywania funkcji jest ujednolicona, pierwsza litera odpowiada formatowi argumentów wejściowych, druga formatowi wyniku następnie występuje sufiks ply(). Litera a oznacza macierz (array), d ramkę danych, l listę, _ nic (funkcji, które nie zwracają wyniku używa się dla ich ,,efektów ubocznych”, np rysowania wykresów). Prefiks r oznacza powtórzenie pewnej operacji n-krotnie na całym zbierze danych, a prefiks m oznacza powtórzenie tej samej funkcji dla zadanego zbioru różnych argumentów wejściowych.
W przypadku pierwszy trzech wierszy tabeli [plyrTabela] funkcje prezentowane różnią się typem wejścia i wyjścia. Mechanizm działania jest podobny. Poniżej przedstawimy tylko funkcje ddply() i llply(). Argumenty tych funkcji są następujące: zbiór danych wejściowych, wskazanie zmiennych grupujących, wskazanie funkcji do zastosowania na każdym podzbiorze, dodatkowe argumenty dla tej funkcji. Zmienne grupujące można wskazać na różne sposoby: podając wektor nazw, formułę lub używając funkcji .(), która nic nie robi ale umożliwia podanie jako argumentów nazw kolumn ze zbioru danych.
W przykładzie poniżej zbiór danych wypadkiRoztopiona podzielimy ze względu na podzbiory określone zmienną geo (czyli podzbiory danych to dane dla różnych krajów). Następnie na każdym podzbiorze wykonamy funkcję lm(), która dla każdego podzbioru wyznaczy współczynniki regresji liniowej – zależności liczby wypadków od roku. Zarówno wejściem jak i wyjściem funkcji jest ramka danych.
wynik <- ddply(wypadkiRoztopiona, .(geo), function(x) lm(value ~ as.numeric(time), data=x)$coef) head(wynik) ## geo (Intercept) as.numeric(time) ## 1 AT 186.6993 -6.097007 ## 2 BE 180.1503 -4.793602 ## 3 BG 142.1307 -1.054696 ## 4 CY 212.8039 -5.646027 ## 5 CZ 158.0719 -2.223942 ## 6 DE 139.9216 -4.868937 |
Często używając funkcji z rodziny **ply() wykorzystuje się funkcje summarize() i transform(). Obie pozwalają na uniknięcie potrzeby tworzenia funkcji anonimowych. Funkcja summarize() jako wynik zwraca wskazane podsumowania każdego podzbioru danych (w poniższym przypadku będą to współczynniki modelu liniowego i średnia liczba wypadków). Funkcja transform() jako wynik zwraca przekształcony zbiór danych z dodanymi nowymi kolumnami.
wynik <- ddply(wypadkiRoztopiona, .(geo), summarize, zmiana = lm(value ~ as.numeric(time))$coef[2], pprzeciecia = lm(value ~ as.numeric(time))$coef[1], srednia = mean(value, na.rm=T)) head(wynik) ## geo zmiana pprzeciecia srednia ## 1 AT -6.097007 186.6993 128.77778 ## 2 BE -4.793602 180.1503 134.61111 ## 3 BG -1.054696 142.1307 132.11111 ## 4 CY -5.646027 212.8039 159.16667 ## 5 CZ -2.223942 158.0719 136.94444 ## 6 DE -4.868937 139.9216 93.66667 |
Poniżej zaprezentujemy przykład użycia funkcji llply(). Oczekuje ona jako argumentu wejściowego listy. Aby taką listę otrzymać posłużymy się funkcją split().
lista <- split(wypadkiRoztopiona, wypadkiRoztopiona$geo) head(lista[[1]]) ## geo value time ## X1991 AT 201 1991 ## X1992 AT 180 1992 ## X1993 AT 163 1993 ## X1994 AT 169 1994 ## X1995 AT 152 1995 ## X1996 AT 129 1996 |
Użycie funkcji llply() jest identyczne z funkcją ddply(), z tą różnicą, że wejście i wyjście jest teraz listą.
wynik <- llply(lista, function(x) lm(value ~ as.numeric(time), data=x)$coef) wynik[[1]] ## (Intercept) as.numeric(time) ## 186.699346 -6.097007 |
Ciekawą funkcją, którą omówimy poniżej jest r*ply(). Powtarza ona wykonanie pewnej operacji zadaną liczbę razy, przez co działa podobnie co funkcja replicate() ale mamy większą kontrolę nad typem wyniku.
wierszy <- nrow(wypadkiRoztopiona) wynik <- rdply(100, lm(value ~ as.numeric(time), data=wypadkiRoztopiona[sample(wierszy, replace=TRUE),])$coef) head(wynik, 3) ## .n (Intercept) as.numeric(time) ## 1 1 169.8856 -4.315885 ## 2 2 178.8349 -4.663816 ## 3 3 171.0432 -4.247858 summary(wynik) ## .n (Intercept) as.numeric(time) ## Min. : 1.00 Min. :162.5 Min. :-5.816 ## 1st Qu.: 25.75 1st Qu.:168.6 1st Qu.:-4.736 ## Median : 50.50 Median :172.5 Median :-4.437 ## Mean : 50.50 Mean :173.0 Mean :-4.465 ## 3rd Qu.: 75.25 3rd Qu.:176.2 3rd Qu.:-4.063 ## Max. :100.00 Max. :187.9 Max. :-3.567 |
Fajnie, a kiedy będzie na rynku?
Można by też uwzględnić pakiet dplyr (https://github.com/hadley/dplyr):
„Currently, it’s not a good idea to have both dplyr and plyr loaded. This is just a short-term problem: in the long-term, I’ll move the matching functions from plyr into dplyr, and add a dplyr dependency to plyr.”
Dzięki, pomysł przyszłościowy. Na dzień dzisiejszy pakiet dplyr jeszcze nie jest dopracowany (może mu się przydarzyć to co pakietowi reshape). Ale jak będzie dostępny w stabilnej wersji to go opiszemy. Na useR HW przedstawiał zyski w czasie przetwarzania i faktycznie dplyr jest znacznie szybszy niż plyr, choć nieznacznie wolniejszy od data.table.
Dla dużych danych bardziej przyszłościowe jednak wydaje mi się bliższa integracja z bazą danych czy hadoopem. Na razie RAM jest zbyt dużym ograniczeniem.
A co do dostępności, to mam nadzieję, że będzie do kupienia od przełomu sierpnia/września. Kończę właśnie skład. Jeszcze czyszczenie z literówek i gotowe.
Możesz zdradzić co jeszcze dodałeś? 🙂
Sporo drobiazgów, ale z dużych rzeczy to opis pakietów knitr, shiny i slidify.
W czwartek będzie wpis o shiny [ta biblioteka jest fantastyczna], w przyszły wtorek o zmianach w typografi [tych jest dużo] a w przyszły czwartek przedmowa do trzeciego wydania.
A jest szansa żeby książka pojawiła się w wersji .mobi?
Jakiś czas temu byłem dużym zwolennikiem wersji mobilnych, ale ograniczała mnie technologia.
Teraz mam wrażenie, że to dobra platforma dla książek liniowych, które można czytać akapit po akapicie, w których nie jest zbyt wielu tabel czy rysunków, skorowidz też nie jest zbyt potrzebny. [jeżeli się mylę to popraw mnie].
,,Przewodnik” ma za to wiele komentarzy na maginesie, odsyłaczy do innych rozdziałów, wtrąceń, uwag. Przerobienie go na taką liniową postać nie byłoby zbyt proste.
Krótka odpowiedź na Twoje pytanie jest więc: może w przyszłości, ale nieprędko.
Kiedy można będzie kupić nowe, trzecie wydanie „Przewodnika…”?
Trudne pytanie. W lipcu myślałem, że będzie to możliwe w sierpniu. Dzisiaj myślę, że będzie to możliwe w grudniu. Napiszę na blogu gdy trzecie wydanie trafi do księgarń.