Pakiet plyr, czyli co nowego w ,,Przewodniku po pakiecie R” wydanie trzecie.

This post was kindly contributed by SmarterPoland » R - go there to comment and to read the full post.

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/paperhttp://www.jstatsoft.org/v40/i01/paper.

Zacznijmy od przedstawienia pakietu reshape2. W tym pakiecie godne uwagi są trzy funkcje: reshape2::melt(), reshape2::acast()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 &#39;http://epp.eurostat.ec.europa.eu/NavTree_prod/everybody/BulkDownloadListing?sort=1'
 file=data%2Ftsdtr420.tsv.gz&#39;
## Content type &#39;application/x-gzip&#39; 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&#39;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()

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