Rozkłady liczebności
Moja dzisiejsza notatka stanowi dopełnienie dwóch poprzednich, w których przypominałem zasady tworzenia rozkładów liczebności i częstości i pokazywałem, jak pisać własne funkcje. Pisanie własnych funkcji ma sens wówczas, gdy jakiś rodzaj analizy wykonujemy bardzo często. Dla przykładu, gdybyśmy chcieli obliczyć wartość indeksu rozbieżności oraz odsetek na głównej przekątnej dla jednej tylko tabeli, być może nie byłoby potrzeby pisania osobnej funkcji. Jeśli jednak często powtarzamy takie obliczenia lub jeśli dysponujemy danymi z bardzo dużej liczby krajów i dla każdego z krajów chcemy wyznaczyć wartosci tych statystyk, wówczas wcześniejsze napisanie odpowiedniej funkcji może bardzo uprościć nam pracę. W tej notatce chcę pokazać — w jaki sposób.
Na początek przypominam funkcję mdi(), której „działanie” opisywałem w poprzedniej notatce. Przypomnę tylko, że funkcja ta (a) stosuje się do macierzy kwadratowych o nieujemnych liczebnościach, (b) oblicza wartości odsetka na głównej przekątnej oraz indeksu rozbieżności między rozkładami brzegowymi, (c) pozwala określić, za pomocą argumentu which, którą z tych wartości funkcja ma zwrócić.
mdi <- function(mat, which = "obie") {
if (which %in% c("główna przekątna", "indeks rozbieżności",
"obie")) {
if (nrow(mat) == ncol(mat)) {
md <- sum(diag(mat))/sum(mat)
di <- sum(abs(rowSums(mat) - colSums(mat)))/(2 * sum(mat))
out <- c(`główna przekątna` = md, `indeks rozbieżności` = di)
if (which == "obie")
out else out[which]
} else cat("Uwaga! Tylko kwadratowe macierze!")
} else cat("Uwaga! argument which wymaga jednej z trzech wartości: 'obie', 'indeks rozbieżności', 'główna przekątna'.")
}We wcześniejszej notatce pokazywałem też, jak za pomocą funkcji xtabs() sporządzić rozkład małżeństw według wykształcenia dla każdego z krajów uczestniczących w piątej rundzie Europejskiego Sondażu Społecznego w 2010 r. Polecenie pokazane w okienku poniżej tworzy tablicę trójwymiarową, której wiersze odpowiadają kategorii wykształcenia męża, kolumny — kategorii wykształcenia żony, zaś trzeci wymiar wskazuje kraj:
tab <- xtabs(pspwght ~ eisced + eiscedp + cntry, data = p, subset = gndr ==
"male" & icpart2 == "yes") + xtabs(pspwght ~ eiscedp + eisced + cntry,
data = p, subset = gndr == "female" & icpart2 == "yes")Przypuśćmy, że chcemy teraz obliczyć odsetki na głównej przekątnej oraz indeksy rozbieżności dla każdego z tych krajów. Gdybyśmy nie dysponowali funkcją mdi(), zadanie to wymagałoby od nas wielokrotnego stosowania tych samych działań do kolejnych „kawałków” tablicy, co jest nie tylko czasochłonne i kłopotliwe, ale naraża nas na częśte błędy. Dla przykładu, pierwszym (w uporządkowaniu alfabetycznym) krajem, jaki obejmuje nasza tabela, jest Albania. Aby obliczyć odsetek na głównej przekątnej dla Albanii, musielibyśmy wpisać:
sum(diag(tab[, , 1]))/sum(tab[, , 1])## [1] 0.64Aby obliczyć analogiczną wielkość dla Belgii, która jest kolejnym krajem w naszej tabeli, musielibyśmy wpisać:
sum(diag(tab[, , 2]))/sum(tab[, , 2])## [1] 0.55I tak dalej, i tak dalej
O pomyłkę łatwo, bo za każdym razem musimy dwukrotnie zmienić wartość indeksu w nawiasie kwadratowym; jeśli o tym zapomnimy, otrzymamy błędny wynik. Dzięki funkcji mdi() zadanie jest łatwiejsze, ale w dalszym ciągu musimy pamiętać, że przed obliczeniem stosowanej wartości dla kolejnego kraju, należy zmienić indeks w nawiasie kwdratowym:
mdi(tab[, , 1], which = "główna przekątna")## główna przekątna
## 0.64mdi(tab[, , 2], which = "główna przekątna")## główna przekątna
## 0.55Inny problem polega na tym, że postępując w ten sposób otrzymamy 29 (bo tyle krajów wzięło udział w 5. rundzie ESS) wartości liczbowych i wykorzystanie ich w dalszych analizach wymagałoby od nas dodatkowych zabiegów. Na przykład, zapisania każdej wartości pod osobną nazwą albo połączenia wszystkich wartości w jeden wektor liczbowy, itd., itp., co znowu naraża nas na błąd.
Jakimś wyjściem z sytuacji jest posłużenie się pętlą for(); w tym celu tworzymy nowy pusty wektor — nazwijmy go md — o długości odpowiadającej liczbie krajów uczestniczących w 5. rundzie ESS, a następnie za pomocą pętli for() podstawiamy wartości odsetka na głównej przekątnej dla kolejnych krajów pod kolejne elementy tego wektora:
md <- vector(mode = "numeric", length = 29)
for (i in 1:29) md[i] <- mdi(tab[, , i], which = "główna przekątna")
md## [1] 0.64 0.55 0.68 0.58 0.61 0.68 0.54 0.56 0.59 0.62 0.58 0.57 0.53 0.67
## [15] 0.59 0.64 0.55 0.61 0.60 0.54 0.58 0.59 0.72 0.67 0.57 0.61 0.70 0.65
## [29] 0.41W analogiczny sposób postępujemy, aby otrzymać zbiór wartości indeksu rozbieżności:
di <- vector(mode = "numeric", length = 29)
for (i in 1:29) di[i] <- mdi(tab[, , i], which = "indeks rozbieżności")
di## [1] 0.081 0.016 0.097 0.175 0.029 0.044 0.168 0.085 0.115 0.042 0.074
## [12] 0.088 0.055 0.046 0.067 0.079 0.075 0.033 0.077 0.060 0.018 0.135
## [23] 0.044 0.055 0.084 0.116 0.066 0.069 0.287Rozwiązanie to daje poprawny wynik, ale mimo wszystko nie jest szczególnie wygodne. Po pierwsze, wymaga od nas utworzenia dwóch osobnych obiektów. Po drugie, wartości każdego z obiektów są nieopisane, jeśli więc chcielibyśmy sprawdzić, ile wynosi wartość indeksu rozbieżności dla, na przykład, Norwegii, musielibyśmy najpierw sprawdzić, na którym miejscu w uporządkowaniu alfabetycznym krajów uczestniczących w 5. rundzie ESS znajduje się Norwegia.
Zaproponujmy zatem inne rozwiązanie, podobne do opisanego wyżej, bo również wykorzystujące pętlę, ale różniące się w ten sposób, że zamiast dwóch osobnych wektorów tworzy jedną macierz o dwóch kolumnach i 29 wierszach; dodatkowo, opatrujemy macierz etykietami, które sprawiają, że jej wartość jest czytelniejsza:
md_di <- matrix(0, ncol = 2, nrow = 29, dimnames = list(Kraj = dimnames(tab)[[3]],
Statystyka = c("główna przekątna", "indeks rozbieżności")))
for (i in 1:29) md_di[i, ] <- mdi(tab[, , i], which = "obie")
md_di## Statystyka
## Kraj główna przekątna indeks rozbieżności
## AL 0.64 0.081
## BE 0.55 0.016
## BG 0.68 0.097
## CH 0.58 0.175
## CY 0.61 0.029
## CZ 0.68 0.044
## DE 0.54 0.168
## DK 0.56 0.085
## EE 0.59 0.115
## ES 0.62 0.042
## FI 0.58 0.074
## FR 0.57 0.088
## GB 0.53 0.055
## HU 0.67 0.046
## IE 0.59 0.067
## IL 0.64 0.079
## IS 0.55 0.075
## IT 0.61 0.033
## LT 0.60 0.077
## NL 0.54 0.060
## NO 0.58 0.018
## PL 0.59 0.135
## PT 0.72 0.044
## RU 0.67 0.055
## SE 0.57 0.084
## SI 0.61 0.116
## SK 0.70 0.066
## UA 0.65 0.069
## XK 0.41 0.287Jak widzimy, argument which funkcji mdi() ma wartość 'obie', co oznacza, iż prosimy funkcję, aby zwróciła wartości obu statystyk: odsetka na głównej przekątnej i indeksu rozbieżności, a następnie podstawiamy uzyskane wartości pod kolejne wiersze macierzy md_di, które odpowiadają kolejnym krajom.
Możemy jednak zaproponować jeszcze prostsze rozwiązanie. Zauważmy bowiem, że zadanie, jakie sobie stawiamy, polega w istocie na:
podziale tablicy
tabna mniejsze „kawałki”, z których każdy reprezentuje rozkład łączny kategorii wykształcenia małżonków w poszczególnym kraju uczestniczących w piątej rundzie ESSzaaplikowaniu funkcji
mdi()do każdego kawałka z osobnapołączeniu wartości zwracanych przez funkcję
mdi()w nowy obiekt.
Zamiast korzystać z pętli, jak w przykładach powyżej, możemy do wykonania tego zadania użyć funkcję apply(), którą zaprojektowano właśnie po to, by wykonywać zadania o opisanej strukturze. Funkcja ta ma następującą składnię:
apply(X = tab, MARGIN = 3, FUN = mdi)Argument X funkcji apply oznacza macierz lub tablicę wielowymiarową (obiekt typu array) lub „ramkę danych” (tzn. obiekt klasy data.frame), na której chcemy wykonać nasze obliczenia. Argument MARGIN odnosi się z kolei do tego, w jaki sposób tablica ma być podzielona na mniejsze elementy; stosują się tutaj te same zasady, które obowiązują w przypadku funkcji prop.table(): w przypadku macierzy MARGIN = 1 aplikuje funkcję do wierszy macierzy, MARGIN = 2 — do jej kolumn, itd. W przykładzie widocznym w okienku powyżej mamy MARGIN = 3, bo tablica, którą utworzyliśmy wcześniej ma 3 wymiary i jej trzeci wymiar odnosi się do kraju. Wreszcie, FUN oznacza funkcję, którą chcemy użyć do obiektu określonego przez X. Powyższe polecenie zwraca:
## cntry
## AL BE BG CH CY CZ DE DK EE
## główna przekątna 0.638 0.548 0.684 0.58 0.607 0.677 0.54 0.560 0.59
## indeks rozbieżności 0.081 0.016 0.097 0.18 0.029 0.044 0.17 0.085 0.11
## cntry
## ES FI FR GB HU IE IL IS
## główna przekątna 0.624 0.577 0.571 0.525 0.672 0.592 0.639 0.548
## indeks rozbieżności 0.042 0.074 0.088 0.055 0.046 0.067 0.079 0.075
## cntry
## IT LT NL NO PL PT RU SE SI
## główna przekątna 0.609 0.601 0.54 0.580 0.59 0.715 0.671 0.568 0.61
## indeks rozbieżności 0.033 0.077 0.06 0.018 0.13 0.044 0.055 0.084 0.12
## cntry
## SK UA XK
## główna przekątna 0.705 0.652 0.41
## indeks rozbieżności 0.066 0.069 0.29Otrzymany obiekt jest macierzą, której wiersze i kolumny są odpowiednio opisane, tak, iż jej zawartość jest czytelna. Jedyny kłopot polega na tym, że nazwy krajów znajdują się w kolumnach, a wartości poszukiwanych statystyk — w wierszach, choć wolelibyśmy, aby było odwrotnie. Możemy, oczywiście, dokonać transpozycji uzyskanej macierzy za pomocą funkcji t():
t(apply(X = tab, MARGIN = 3, FUN = mdi))##
## cntry główna przekątna indeks rozbieżności
## AL 0.64 0.081
## BE 0.55 0.016
## BG 0.68 0.097
## CH 0.58 0.175
## CY 0.61 0.029
## CZ 0.68 0.044
## DE 0.54 0.168
## DK 0.56 0.085
## EE 0.59 0.115
## ES 0.62 0.042
## FI 0.58 0.074
## FR 0.57 0.088
## GB 0.53 0.055
## HU 0.67 0.046
## IE 0.59 0.067
## IL 0.64 0.079
## IS 0.55 0.075
## IT 0.61 0.033
## LT 0.60 0.077
## NL 0.54 0.060
## NO 0.58 0.018
## PL 0.59 0.135
## PT 0.72 0.044
## RU 0.67 0.055
## SE 0.57 0.084
## SI 0.61 0.116
## SK 0.70 0.066
## UA 0.65 0.069
## XK 0.41 0.287Istnieje jednak jeszcze prostsze rozwiązanie, które polega na posłużeniu się funkcją adply() z pakietu plyr. Cały pakiet stanowi uogólnienie podejścia do analizy danych opartego na podziale analizowanego obiektu na mniejsze części, zaaplikowaniu odpowiedniej funkcji do każdego kawałka i połączeniu otrzymanych rezultatów w nowy obiekt (split - apply - combine). Uogólnieniem, ponieważ różne funkcje służące realizacji tego celu są wbudowane w R, np. opisana wyżej funkcją apply(), ale też inne: tapply(), sapply() czy aggregate(). Funkcje te różnią się typem obiektów, do jakich można je zastosować — niektóre stosuje się do macierzy lub tablic, inne do list lub ramek danych — oraz rodzajem obiektu, który funkcja zwraca. Pakiet plyr zawiera całą rodzinę funkcji, które łączą ostatnie trzy litery w nazwie — ply — zaś dwie pierwsze służą do określenia, jakiego typu jest obiekt wsadowy i jakiej klasy ma być obiekt „wyjściowy”. I tak, na przykład, funkcja aaply() wymaga tablicy (array) i zwraca tablicę, funkcja daply() wymaga ramki danych (data.frame), lecz zwraca tablicę, funkcja ldply() wymaga listy i zwraca ramkę danych, itd., itp. Więcej szczegółów na temat poszczególnych funkcji w pakiecie plyr, a także, hm, filozofii, na jakiej pakiet ten się opiera, można znaleźć w tym artykule.
Wspomniana wcześniej funkcja adply() wymaga, zgodnie z tą konwencją, tablicy i zwraca ramkę danych. Mówiąc dokładniej, funkcja ta dzieli tablicę wejściową na mniejsze kawałki w sposób określony przez użytkownika, do każdego kawałka aplikuje funkcję wskazaną przez użytkownika, a następnie łączy wyniki obliczeń w ramkę danych. Dla danych z naszego przykładu mamy więc:
library(plyr)
adply(.data = tab, .margin = 3, .fun = mdi)## cntry główna przekątna indeks rozbieżności
## 1 AL 0.64 0.081
## 2 BE 0.55 0.016
## 3 BG 0.68 0.097
## 4 CH 0.58 0.175
## 5 CY 0.61 0.029
## 6 CZ 0.68 0.044
## 7 DE 0.54 0.168
## 8 DK 0.56 0.085
## 9 EE 0.59 0.115
## 10 ES 0.62 0.042
## 11 FI 0.58 0.074
## 12 FR 0.57 0.088
## 13 GB 0.53 0.055
## 14 HU 0.67 0.046
## 15 IE 0.59 0.067
## 16 IL 0.64 0.079
## 17 IS 0.55 0.075
## 18 IT 0.61 0.033
## 19 LT 0.60 0.077
## 20 NL 0.54 0.060
## 21 NO 0.58 0.018
## 22 PL 0.59 0.135
## 23 PT 0.72 0.044
## 24 RU 0.67 0.055
## 25 SE 0.57 0.084
## 26 SI 0.61 0.116
## 27 SK 0.70 0.066
## 28 UA 0.65 0.069
## 29 XK 0.41 0.287W efekcie naszych działań dostajemy nowy zbiór danych, który możemy wykorzystać w kolejnych analizach. Na przykład:
tab_2 <- adply(.data = tab, .margin = 3, .fun = mdi)
plot(tab_2$"główna przekątna" ~ tab_2$"indeks rozbieżności",
las = 1, pch = 19, col = "navyblue", main = "Regresja liniowa",
xlab = "Indeks rozbieżności", ylab = "Odsetek na głównej przekątnej")
abline(lm(tab_2$"główna przekątna" ~ tab_2$"indeks rozbieżności"),
col = "firebrick", lwd = 2)
Rysunek powyżej pokazuje, jak zmieniają się wartości jednej zmiennej wraz ze zmianą wartości drugiej z nich. Na diagram rozrzutu nałożono linię regresji. Widzimy, iż wartość odsetka na głównej przekątnej maleje wraz ze wzrostem indeksu rozbieżności. To, oczywiście, tylko bardzo prosty przykład; można jednak zaproponować przykłady znacznie bardziej złożone, w których liczba jednostek analizy jest daleko większa niż w naszym przykładzie, a i liczba zmiennych, które uzyskuje się w wyniku zastosowanych przekształceń jest dużo większa.
Zdaję sobie sprawę, że pointa notki dzisiejszej i dwóch poprzednich może niknąć w całym tym moim gadulstwie, chcę zatem podkreślić, że sedno wszystkich podawanych przeze mnie przykładów i wyjaśnień sprowadza się do tego, że od zbioru danych surowych do wykresu pokazanego wyżej potrzeba tylko kilku kroków:
Po pierwsze, tworzymy tablicę trójwymiarową zawierającą rozkłady łączne kategorii wykształcenia małżonków w krajach uczestniczących w ESS:
tab <- xtabs(pspwght ~ eisced + eiscedp + cntry, data = p, subset = gndr ==
"male" & icpart2 == "yes") + xtabs(pspwght ~ eiscedp + eisced + cntry, data = p,
subset = gndr == "female" & icpart2 == "yes")Po drugie, tworzymy funkcję, która ma zwracać wartości poszukiwanych przez nas statystyk:
mdi <- function(mat, which = "obie") {
if (which %in% c("główna przekątna", "indeks rozbieżności",
"obie")) {
if (nrow(mat) == ncol(mat)) {
md <- sum(diag(mat))/sum(mat)
di <- sum(abs(rowSums(mat) - colSums(mat)))/(2 *
sum(mat))
out <- c(`główna przekątna` = md, `indeks rozbieżności` = di)
if (which == "obie")
out else out[which]
} else cat("Uwaga! Tylko kwadratowe macierze!")
} else cat("Uwaga! argument which wymaga jednej z trzech wartości: 'obie', 'indeks rozbieżności', 'główna przekątna'.")
}Po trzecie, aplikujemy tę funkcję do każdej z tabel cząstkowych, reprezentujących poszczególne kraje:
adply(.data = tab, .margin = 3, .fun = mdi)I tyle. I to już wszystko. Ostatnie polecenie daje nam ramkę danych, którą wykorzystamy do wykresu.