Chciałbym dzisiaj na kilku prostych przykładach pokazać możliwości R w zakresie analizy danych pochodzących z międzynarodowych badań społecznych czy badań edukacyjnych. W pierwszym przypadku możemy być zainteresowani poziomem nierówności zarobków w różnych krajach, szacowanym na podstawie odpowiedzi na pytanie o przeciętne miesięczne zarobki udzielanych przez respondentów biorących udział w badaniu sondażowym. W badaniach edukacyjnych często chcemy ocenić różnice między szkołami pod względem średniego wyniku uczniów w teście umiejętności. Zestawienie to możemy łatwo uzyskać odpowiednio modyfikując zamieszczone w tej notce polecenia.
W pewnym badaniu socjologicznym respondentów zapytano o to, ile według nich zarabiają przedstawiciele kilku zawodów, oraz o to, ile powinni oni zarabiać. Zebrane dane zapisno w zbiorze danych o następującej postaci:
## respid occupation act.earn just.earn
## 1 1 CEO 11000 10000
## 2 1 GP in a public hospital 2500 3500
## 3 1 attorney 7500 5000
## 4 1 bricklayer 1000 2000
## 5 1 city bus driver 1500 2000
## 6 1 small-shop owner 2000 2500
## 7 2 bank clerk 2000 3000
## 8 2 GP in a public hospital 3000 6000
## 9 2 farm worker 2000 4000
## 10 2 bricklayer 4000 3500
## 11 2 city bus driver 3000 5000
## 12 2 small-shop owner 2000 3000
Jak widzimy, w zbiorze znajdują się 4 zmienne:
respid, która oznacza identyfikator respondentaoccupation, zawierająca nazwy zawodów, których zarobki respondenci ocenialiact.earn, zawierająca oceny rzeczywistych zarobków w ocenianaych zawodachjust.earn, na którą składają się oceny zarobków sprawiedliwych w tych zawodach.De facto, w badaniu mierzono również inne zmienne, ponieważ jednak nie będziemy z nich korzystać w omawianych niżej przykładach, ograniczymy się tylko do zaprezentowanych 4 zmiennych.
W badaniu udział wzięło łącznie \( n=120 \) respondentów i każdy z nich oceniał zarobki w \( m=6 \) zawodach, co daje łącznie 720 obserwacji. Jeden z respondentów nie udzielił odpowiedzi na pytanie o ocenę zarobków rzeczywistych; respondentowi temu przypisano wartość mediany ocen wszystkich pozostałych respondentów. Dodatkowo, niektórzy respondenci podawali nie pojedyncze kwoty, lecz przedziały, w jakich mieszczą się według nich zarobki rzeczywiste i sprawiedliwe w ocenianych zawodach. Odpowiedzi takie zastępowano kwotą stanowiącą środek przedziału.
Naszym celem jest oszacowanie poziomu nierówności postrzeganej, tj. nierówności ocen zarobków rzeczywistych, oraz nierówności postulowanej, tj. nierówności ocen zarobków sprawiedliwych, dla każdego respondenta. Wymaga to od nas:
act.earn i just.earn funkcji zwracającej poziom nierównościDo oceny poziomu nierówności zarobków wykorzystamy współczynnik Giniego. W R nie ma wbudowanej funkcji zwracającej wartość współczynnika Giniego dla zadanego zbioru wartości, można ją jednak znaleźć w pakiecie ineq, można też napisać ją samemu. Skorzystamy z tej ostatniej opcji. Posłużymy się w tym celu następującym wzorem: \[ G(X)=\frac{N+1}{N}-\frac{2}{E(X)N^2}\sum_{j=1}^N(N+1-j)x_j \] We wzorze tym zakłada się, że wartości zmiennej \( X \) są uporządkowane niemalejąco: \( x_1 \leq x_2 \leq \cdots \leq x_j \leq \cdots \leq x_N \). Funkcja, która oblicza wartość współczynnika Giniego według tego wzoru, ma następującą postać:
gini <- function(x, na.rm=TRUE) {
if (na.rm) x <- x[ !is.na(x) ] else if ( any(is.na(x)) ) return(NA)
x <- sort(x)
N <- length(x)
j <- order(x)
(N+1)/N - 2/( mean(x)*N^2 )*sum( (N+1-j)*x )
}
Jak widzimy, funkcja ta wymaga dwóch argumentów: x, czyli wektora wartości liczbowych, dla których chcemy obliczyć wartość współczynnika Giniego, oraz na.rm, czyli wartości logicznej, która określa sposób postępowania w przypadku występowania w zbiorze x braków danych. Jeśli na.rm==TRUE, wówczas braki danych, jeśli występują, zostają wykluczone ze zbioru x i pominięte przy wykonywaniu obliczeń. W przeciwnym przypadku, tzn. gdy na.rm==FALSE, zaś x zawiera braki danych, program zwróci NA. Domyślnie, wartość argumentu na.rm jest ustawiona na TRUE.
Po zdefiniowaniu funkcji gini możemy przystąpić do rozwiązywania postawionego wyżej problemu.
Pierwsze rozwiązanie tego problemu polega na
act.earn w macierz o wymarach \( 120 \times 6 \), której wiersze ,,odpowiadają'' respondentom, zaś kolumny — zawodom, których zarobki są oceniane, a następniegini(x, na.rm)just.earn.aux.act <- matrix(p$act.earn, ncol = 6, byrow = TRUE)
act.gini <- apply(X = aux.act, MARGIN = 1, FUN = gini)
summary(act.gini)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0409 0.1620 0.2770 0.2870 0.3810 0.6970
aux.just <- matrix(p$just.earn, ncol = 6, byrow = TRUE)
just.gini <- apply(X = aux.just, MARGIN = 1, FUN = gini)
summary(just.gini)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0291 0.1360 0.1950 0.2130 0.2610 0.6460
Zaproponowane rozwiązanie jest w gruncie rzeczy proste, ale wymaga od nas utworzenia dwóch dodatkowych, pomocniczych zbiorów danych, aux.act i aux.just, z których po zakończeniu obliczeń, nie będziemy już korzystać. Ponadto, uzyskujemy w efekcie dwie osobne zmienne — act.gini i just.gini — a nie zbiór danych, o jaki nam chodziło.
Poszukajmy zatem innego rozwiązania.
Drugie rozwiązanie polega na podziale, za pomocą funkcji split wyjściowego zbioru danych na 120 ,,kawałków'', z których każdy zawiera tylko tę część zbioru wyjściowego, w której znajdują się odpowiedzi poszczególnego respondenta:
aux <- split(p[, c("act.earn", "just.earn")], p$respid)
Za pomocą powyższego polecenia tworzymy listę składającą się ze 120 elementów, przy czym każdy element zawiera podzbiór wyjściowego zbioru danych, w którym znajdują się tylko dane dotyczące poszczególnego respondenta. Zauważmy także, że utworzenia tej listy wykorzystaliśmy jedynie dwie zmienne — oceny zarobków rzeczywistych i sprawiedliwych — a nie wszystkie cztery zmienne zawarte w oryginalnym zbiorze.
Zauważmy teraz, że pojedynczy element utworzonej listy możemy potraktować jak macierz o wymiarach \( 6\times 2 \), której kolumny zawierają oceny zarobków, odpowiednio, rzeczywistych i sprawiedliwych podane przez poszczególnego respondenta. Aplikując funkcję gini do kolumn tej macierzy otrzymamy więc poziom nierówności postrzeganej i postulowanej przez tego respondenta. Dla przykładu, w przypadku pierwszego respondenta otrzymujemy:
apply(X = aux[[1]], MARGIN = 2, FUN = gini)
## act.earn just.earn
## 0.4477 0.3333
Funkcja apply wymaga określenia wartości trzech argumentów:
X to nazwa obiektu, który chcemy poodać przekształceniu, obiektem tym najczęściej będzie macierz lub zbiór danychFUN to funkcja, którą chcemy zastosować. MARGIN określa ,,kierunek'' wykonywania obliczeń; gdy MARGIN==1, funkcja jest aplikowana do wierszy, gdy zaś MARGIN==2, jest ona aplikowana do kolumn.Obok nazwy funkcji konieczne jest czasami podanie dodatkowych informacji lub argumentów, wymaganych przez funkcję określoną przez FUN. Dla przykładu, gdybyśmy chcieli obliczyć dla danego respondenta wartość pierwszego i trzeciego kwartyla ocen zarobków postrzeganych i postulowanych, musielibyśmy użyć następującego polecenia:
apply(X = aux[[1]], MARGIN = 2, FUN = quantile, probs = c(0.25, 0.75))
## act.earn just.earn
## 25% 1625 2125
## 75% 6250 4625
W tym przykładzie, argument probs jest argumentem funkcji quantile, określającym, które kwantyle chcemy obliczyć, nie zaś argumentem funkcji apply.
Wróćmy teraz do naszego podstawowego problemu, a mianowicie obliczenia poziomu nierówności w ocenach zarobków postrzeganych i postulowanych. Zależy nam jednak na obliczeniu tych wielkości dla wszystkich respondentów w próbie, funkcję apply musimy więc zastosować do wszystkich elementów listy aux. Najwygodniej w tym celu skorzystać z funkcji sapply, która aplikuje zadaną funkcję do wszystkich elementów listy:
ds <- sapply(X = aux, FUN = function(x) apply(X = x, MARGIN = 2, FUN = gini))
Uzyskujemy w ten sposób macierz o wymiarach \( 2\times 120 \). Za pomocą funkcji t() dokonamy transpozycji tej macierzy do pożądanej postaci:
ds <- t(ds)
head(ds)
## act.earn just.earn
## 1 0.44771 0.3333
## 2 0.14583 0.1463
## 3 0.09211 0.1887
## 4 0.16667 0.1768
## 5 0.36032 0.2925
## 6 0.17126 0.1389
Rozwiązanie to jest nieco wygodniejsze niż poprzednie, ale też mniej intuicyjne — zwłaszcza w przedostatnim kroku. Spróbujmy zatem poszukać jeszcze innego rozwiązania dla naszego problemu.
Rozwiązanie trzecie korzysta z funkcji dostępnych w pakiecie plyr. Pakiet ten jest zbiorem funkcji ułatwiających przekształcanie danych do postaci dogodnej do analizy, w szczególności, w sytuacji, gdy analiza ma na celu (lub obejmuje) podział wyjściowego zbioru danych na mniejsze elementy, zaaplikowanie danej funkcji do każdego elementu i połączenie wartości zwróconych przez tę funkcję w nowy zbiór.
Pierwsza funkcja, z której będziemy korzystać, to funkcja summarise, która przekształca wyjściowe wartości zmiennych w zadany sposób. Dla przykładu, polecenie:
library(plyr)
summarise(p, act.gini = gini(act.earn), just.gini = gini(just.earn))
zwróci nam nowy obiekt, który składał się będzie z jednego tylko wiersza zawierającego informacje o poziomie nierówności postrzeganej i postulowanej dla całej próby:
## act.gini just.gini
## 1 0.4757 0.3803
Pierwszym argumentem funkcji summarise jest nazwa zbioru danych, w którym znajdują się zmienne, jakie zamierzamy poddać przekształceniu. Kolejne argymentu zawierają wyrażenia określające, w jaki sposób zmienne mają być przekształcone.
W podanym przykładzie uzyskujemy poziom nierówności postrzeganej i postulowanej dla całej próby. Nas jednak interesuje poziom nierówności postrzeganej i postulowanej przez każdego z respodentów osobno. Chcemy więc funkcję summarise zastosować osobno do każdego podzbioru wyjściowego zbioru danych, w którym znajdują się odpowiedzi poszczególnego respondenta. Wykorzystamy w tym celu funkcję ddply.
ds <- ddply(.data = p, .variables = "respid", .fun = summarise, act.gini = gini(act.earn),
just.gini = gini(just.earn))
Funkcja ddply ma trzy obowiązkowe argumenty:
summarise wraz z dodatkowymi argumentami określającymi rodzaj przekształceń, którym chcemy poddać zmienne w wyjściowym zbiorze danych.Po wykonaniu powyższego polecenia uzyskamy zbiór danych o następującej postaci:
## respid act.gini just.gini
## 1 1 0.44771 0.3333
## 2 2 0.14583 0.1463
## 3 3 0.09211 0.1887
## 4 4 0.16667 0.1768
## 5 5 0.36032 0.2925
## 6 6 0.17126 0.1389
Jest to zatem dokładnie ten rezultat, o który nam chodziło.
Na rysunku poniżej przedstawiono diagram rozrzutu dla nierówności postrzeganej i postulowanej:
I na zakończenie jeszcze zestawienie wszystkich trzech proponowanych rozwiązań:
## ROZWIĄZANIE PIERWSZE
aux.act <- matrix(p$act.earn, ncol = 6, byrow = TRUE)
act.gini <- apply(X = aux.act, MARGIN = 1, FUN = gini)
summary(act.gini)
aux.just <- matrix(p$just.earn, ncol = 6, byrow = TRUE)
just.gini <- apply(X = aux.just, MARGIN = 1, FUN = gini)
summary(just.gini)
## ROZWIĄZANIE DRUGIE
aux <- split(p[, c("act.earn", "just.earn")], p$respid)
ds <- sapply(X = aux, FUN = function(x) apply(X = x, MARGIN = 2, FUN = gini))
ds <- t(ds)
## ROZWIĄZANIE TRZECIE
library(plyr)
ds <- ddply(.data = p, .variables = "respid", .fun = summarise, act.gini = gini(act.earn),
just.gini = gini(just.earn))