W poprzedniej notatce opisywałem tworzenie rozkładów częstosci i liczebności dwóch lub więcej zmiennych za pomocą funkcji xtabs(). Zwykle rozkłady częstości czy liczebności jako takie nas nie interesują, interesują nas natomiast pewne ich własności, takie jak siła zależności między kategoriami w wierszach i kolumnach tabeli albo, w przypadku tablic „kwadratowych”, odsetek przypadków na głównej przekątnej. Czasami dla jednej tabeli chcemy określić więcej niż jedną własność. Czasami dysponujemy, tak jak w przykładach opisywanych poprzednio, kilkudziesięcioma tabelami i chcemy wyznaczyć różne własności dla każdej tabeli z osobna. Szukamy zatem prostej metody określania tych własności dla wielu tablic jednocześnie.

W niniejszej notatce postaram się pokazać, jak rozwiązywać tego rodzaju problemy posługując się tymi samymi danymi, co poprzednio, a więc danymi zdającymi sprawę z rozkładu łącznego kategorii wykształcenia małżonków (partnerów) w 29 krajach uczestniczących w 6. rundzie Europejskiego Sondażu Społecznego. Dane są publicznie dostępne i każdy może je pobrać ze strony ESS po zarejestrowaniu się w bazie użytkowników danych; procedura rejestracji jest, jak pisałem, bezbolesna.

Rozkłady łączne kategorii wykształcenia małżonków mają postać macierzy kwadratowych, a więc macierzy o jednakowej liczbie wierszy i kolumn. Macierze takie są specyficzne, ponieważ mają pewne własności, których nie da się określić dla macierzy innego rodzaju. Na przykład, dla macierzy kwadratowych można wyznaczyć: (a) odsetek przypadków na głównej przekątnej, (b) zgodność rozkładów kategorii wierszowych i kolumnowych, (c) wielkość rozbieżności między rozkładami brzegowymi, by podać tylko kilka przykładów. Mniejsza jednak o konkretne przykłady, chodzi mi raczej o nakreślenie metody postępowania, która nie ogranicza się bynajmniej do tablic kwadratowych. Dla przykładu, siłę zależności między kategoriami wierszowymi i kolumnowymi można wyznaczyć dla tablic o różnych wymiarach, nie tylko dla tablic kwadratowych.

Odsetki na głównej przekątnej

Niektóre z opisanych wyżej własności tablic można wyznaczyć za pomocą funkcji dostępnych w R. W wielu przypadkach jest jednak tak, że interesuje nas jakaś szczególna własność, dla której nie ma odpowiedniej funkcji lub istniejące funkcje nie są wystarczające. Rozpatrzmy, dla ilustracji, odsetek obserwacji na głównej przekątnej. Zacznijmy od utworzenia tabeli przedstawiającej rozkład łączny kategorii wykształcenia mężów i żon (lub partnerów i partnerek) w Polsce w 2012 r. Korzystam w tym celu z poleceń omawianych w poprzedniej notatce:

tab_1 <- xtabs(pspwght ~ eisced + eiscedp, data = p, subset = gndr == "male" & icpart2 == 
    "yes" & cntry == "PL") + xtabs(pspwght ~ eiscedp + eisced, data = p, subset = gndr == 
    "female" & icpart2 == "yes" & cntry == "PL")
tab_1
##                  eiscedp
## eisced            primary lower secondary upper secondary post-secondary
##   primary             5.9             9.1             0.0            0.0
##   lower secondary     9.9           322.9           160.4           68.1
##   upper secondary     0.0            68.9           153.9          101.7
##   post-secondary      0.0             6.9            41.9          187.8

Uzyskujemy w efekcie macierz o wymiarach 4 × 4. Dostępna w R funkcja diag() zwraca elementy na głównej przekątnej:

diag(tab_1)
##         primary lower secondary upper secondary  post-secondary 
##             5.9           322.9           153.9           187.8

Mówiąc inaczej, funkcja diag() zwraca informację o liczbie par w poszczególnych komórkach na głównej przekątnej, nam tymczasem chodzi odsetek obserwacji na głównej przekątnej, co w kontekście danych na temat małżeństw według poziomu wykształcenia interpretuje się jako odsetek małżeństw, w których mężowie i żony mają ten sam poziom wykształcenia. Aby uzyskać wielkość, której szukamy, musimy dodać do siebie wartości zwracane przez funkcję diag(), a następnie podzielić przez sumę obserwacji w całej tabeli:

sum(diag(tab_1))/sum(tab_1)
## [1] 0.59

Alternatywnie, możemy najpierw tabelę liczebności zamienić na tabelę proporcji, a następnie zsumować wartości na głównej przekątnej:

sum(diag(prop.table(tab_1)))
## [1] 0.59

Jak widzimy, w obu przypadkach dostajemy dokładnie ten sam rezultat. Zaproponowane rozwiązania, jakkolwiek dają poprawny wynik, są jednak mało wygodne, jeśli chcemy obliczyć odsetki na głównej przekątnej dla wielu tabel, a tak jest w przypadku analizy danych z Europejskiego Sondażu Społecznego czy innych projektów porównawczych. Stosowanie zaproponowanej wyżej formuły do każdej z tabeli z osobna jest czasochłonne i wiąże się z wysokim ryzykiem błędu. Ryzyko to można znacznie zmniejszyć poprzez: (a) utworzenie funkcji, która zwraca odsetek obserwacji na głównej przekątnej oraz (b) skorzystanie z licznych w R rozwiązań w zakresie przekształceń grupowych, tzn. rozwiązań, które umożliwiają aplikowanie jednej funkcji do różnych podzbiorów tego samego zbioru jednocześnie. W dalszej części notatki skupiam się na tym pierwszym problemie. W kolejnej notatce opiszę rozwiązania drugiego z nich.

Tworzenie własnych funkcji

No tak, tworzenie własnych funkcji w R to jest temat rzeka… Niemal wszystko, czym posługujemy się w R, jest funkcją. Posługiwanie się funkcjami oznacza, że na wejściu mamy pewien argument lub zbiór argumentów, które przekształcamy za pomocą jednej bądź szeregu operacji, a w wyniku tych przekształceń uzyskujemy nową wartość lub zbiór wartości. W przypadku funkcji, która dla zadanej macierzy kwadratowej ma zwracać odsetek przypadków na głównej przekątnej, argumentem jest właśnie zadana macierz, przekształcenie polega na zsumowaniu elementów na głównej przekątnej i podzieleniu przez sumę wszystkich elementów macierzy, zaś wynik tych obliczeń to wartość, której szukamy. Nazwijmy poszukiwaną przez nas funkcję m_diag(); aby ją utworzyć, posłużymy się następującym poleceniem:

m_diag <- function(mat) sum(diag(mat))/sum(mat)

W powyższym zapisie, mat oznacza zmienną, w miejsce której podstawiamy nazwię macierzy, dla której chcemy obliczyć odsetek obserwacji na głównej przekątnej. Aplikując tę funkcję do utworzonej wcześniej macierzy tab_1, zawierającej informację o liczebnościach małżeństw wg poziomu wykształcenia w Polsce w 2010 r., otrzymujemy:

m_diag(mat = tab_1)
## [1] 0.59

Jak widzimy, nasza funkcja jest bardzo prosta: wymaga pojedynczego argumentu i zwraca pojedynczą wartość. Możemy jednak tworzyć funkcje o bardziej złożonej postaci. Mówienie o „głównej przekątnej” ma sens tylko w przypadku macierzy kwadratowych. Spróbujmy więc zmodyfikować naszą funkcję w taki sposób, aby przed wykonaniem odpowiednich obliczeń sprawdzała, czy macierz „wejściowa” jest kwadratowa, tzn. liczba wierszy jest równa w jej przypadku liczbie kolumn; jeśli jest, funkcja wykona opisane wyżej obliczenia i zwróci ich wynik, w przeciwnym razie — zwróci komunikat o błędzie:

m_diag <- function(mat) {
    if (nrow(mat) == ncol(mat)) 
        sum(diag(mat))/sum(mat) else cat("Uwaga! Tylko kwadratowe macierze!")
}
m_diag(mat = tab_1)
## [1] 0.59
m_diag(mat = tab_1[-1, ])
## Uwaga! Tylko kwadratowe macierze!

Powyższe polecenia modyfikują naszą funkcję do omawianej postaci, a następnie sprawdzają jej działanie na dwóch macierzach: wykorzystywanej wcześniej przez nas macierzy tab_1 oraz tej samej macierzy, ale bez pierwszego wiersza. Usunięcie pierwszego wiersza sprawia, że macierz nie jest już kwadratowa, nie można więc wyznaczyć dla niej głównej przekątnej i w konsekwencji funkcja zwraca komunikat o błędzie.

Możemy teraz rozbudować naszą funkcję jeszcze bardziej. Zauważmy, że suma (bądź odsetek) elementów na głównej przekątnej to nie jedyna charakterystyka macierzy kwadratowej, która może interesować badacza. Inną taką charakterystyką jest różnica między rozkładami brzegowymi. W przypadku tablic reprezentujących łączny rozkład poziomów wykształcenia mężów i żon różnica ta zdaje sprawę z rozbieżności między strukturami wykształcenia mężów i żon, która „wymusza” zawieranie małżeństw z osobami spoza własnej kategorii wykształcenia. Mówiąc inaczej, w zbiorowościach, w których struktury wykształcenia kobiet i mężczyzn odbiegają od siebie, jakaś część małżeństw musi być między osobami o różnym poziomie wykształcenia — bez względu na indywidualne preferencje co do wykształcenia partnera.

Do określenia tej różnicy wykorzystuje się indeks rozbieżności (index of dissimilarity), który stanowi miarę odległości między dwoma rozkładami. Niech \(A\) oznacza macierz rozkładu łącznego poziomów wykształcenia kobiet i mężczyzn. Przez \(a_{gh}\) oznaczamy element w \(g\)-tym wierszu i \(h\)-tej kolumniej macierzy \(A\), natomiast przez \(a_{g+}\) i \(a_{+h}\) sumę elementów w, odpowiednio, \(g\)-tym wierszu i \(h\)-tej kolumnie. Indeks rozbieżności definiuje się następująco: \[ID = \frac{1}{2}\sum_{g=1}^m \lvert a_{g+} - a_{+g} \rvert,\] gdzie \(m\) oznacza liczbę wierszy (kolumn) macierzy \(A\). Przyjmuję tu, że elementy macierzy \(A\) reprezentują odsetki, a nie liczebności.

Wartość indeksu rozbieżności dla zadanej macierzy kwadratowej możemy wyznaczyć za pomocą następującej funkcji:

diss <- function(mat) {
    if (nrow(mat) == ncol(mat)) 
        sum(abs(rowSums(mat) - colSums(mat)))/(2 * sum(mat)) else cat("Uwaga! Tylko kwadratowe macierze!")
}
diss(mat = tab_1)
## [1] 0.13
diss(mat = tab_1[-1, ])
## Uwaga! Tylko kwadratowe macierze!

Funkcja diss() „zachowuje się” w sposób bardzo zbliżony do utworzonej wcześniej funkcji m_diag(): najpierw sprawdza, czy macierz wejściowa jest kwadratowa; jeśli jest, wówczas funkcja oblicza wartość indeksu rozbieżności według podanej wyżej formuły, w przeciwnym razie zwraca komunikat o błędzie.

Zauważmy jednak, że zamiast dwóch osobnych funkcji, możemy utworzyć jedną, która będzie obliczać wartości obu statystyk; nazwijmy tę funkcję mdi() (od main diagonal oraz dissimilariy index):

mdi <- function(mat) {
    if (nrow(mat) == ncol(mat)) {
        md <- sum(diag(mat))/sum(mat)
        di <- sum(abs(rowSums(mat) - colSums(mat)))/(2 * sum(mat))
        print(c(md, di))
    } else cat("Uwaga! Tylko kwadratowe macierze!")
}
mdi(tab_1)
## [1] 0.59 0.13

Tak jak we wcześniejszych przypadkach, utworzona przez nas funkcja najpierw sprawdza, czy macierz jest kwadratowa; jeśli nie jest, zwraca komunikat o błędzie; jeśli jest, oblicza wartości obu poszukiwanych statystyk, a następnie wyświetla ja na ekranie. Mówiąc dokładniej, jeśli macierz jest kwadratowa, funkcja mdi() podstawia pod md wartość odsetka na głównej przekątnej, a pod di wartość indeksu rozbieżności, po czym obie wielkości łączy w jeden wektor, który wyświetla na ekranie. Takie rozwiązanie jest poprawne, aczkolwiek mało wygodne, ponieważ wymaga od użytkownika pamiętania o tym, która z dwu wartości wyświetlanych na ekranie odpowiada której statystyce. Problem ten możemy obejść modyfikując nieco funkcję mdi():

mdi <- function(mat) {
    if (nrow(mat) == ncol(mat)) {
        cat("Odsetek elementów na głównej przekątnej:\n", sum(diag(mat))/sum(mat), 
            "\nIndeks rozbieżności:\n", sum(abs(rowSums(mat) - colSums(mat)))/(2 * 
                sum(mat)))
    } else cat("Uwaga! Tylko kwadratowe macierze!")
}
mdi(tab_1)
## Odsetek elementów na głównej przekątnej:
##  0.59 
## Indeks rozbieżności:
##  0.13
mdi(tab_1[-1, ])
## Uwaga! Tylko kwadratowe macierze!

Jak widzimy, tym razem wartości zwracane przez funkcję mdi() są dokładnie opisane. Zamiast, jak poprzednio, tworzyć dwa nowe obiekty md oraz di, aby podstawić pod nie wartości obu statystyk, wykorzystaliśmy funkcję cat(), która pozwala na łączenie obiektów tekstowych i wyników obliczeń na ekranie; widoczny w poleceniach w powyższym okienku symbol \n oznacza znak nowej linii, dzięki czemu poszczególne elementy wyświetlają się na ekranie jeden pod drugim.

Wartość, jaką zwraca funkcja, zależy w dużej mierze od tego, jaki użytek zamierzamy uczynić z tej wartości. Bardzo często chcemy wykorzystać efekt bieżących obliczeń do sporządzenia wykresu, tabeli czy po prostu w dalszych analizach, ważne jest więc, aby rezultat zwracany przez funkcję miał łatwą do przetwarzania postać. Rezultat widoczny w okienku powyżej łatwy do przetwarzania nie jest, ma on bowiem postać bloku tekstowego, z którego dopiero trzeba wyodrębnić wartości liczbowe za pomocą dodatkowych narzędzi. Możemy jednak nieco zmodyfikować naszą funkcję, tak, by zwracała wartości liczbowe opatrzone odpowiednimi nazwami:

mdi <- function(mat) {
    if (nrow(mat) == ncol(mat)) {
        md <- sum(diag(mat))/sum(mat)
        di <- sum(abs(rowSums(mat) - colSums(mat)))/(2 * sum(mat))
        c(`główna przekątna` = md, `indeks rozbieżności` = di)
    } else cat("Uwaga! Tylko kwadratowe macierze!")
}
mdi(tab_1)
##    główna przekątna indeks rozbieżności 
##                0.59                0.13

Jeśli do dalszych obliczeń chcemy wybrać tylko wartość indeksu rozbieżności, wystarczy, jeśli wpiszemy:

mdi(tab_1)["indeks rozbieżności"]
## indeks rozbieżności 
##                0.13

Możemy też rozszerzyć naszą funkcję o dodatkowy argument określający, którą z dwu wartości funkcja ma zwrócić; po tych zmianach funkcja miałaby następującą postać:

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") 
                print(out) else print(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'.")
}

W porównaniu z poprzednimi wersjami, wersja w okienku powyżej zawiera dwie istotne zmiany. Po pierwsze, jak widzimy w pierwszym wierszu, funkcja wymaga teraz dwóch argumentów, nie jednego; ten nowy argument został nazwany which i określa, którą z dwóch wartości — odsetek na głównej przekątnej czy indeks rozbieżności — funkcja ma zwrócić. Inaczej jednak niż argument mat, który nie ma wartości domyślnej, argument which taką wartość posiada — jak bowiem widzimy argument ten ma w deklaracji funkcji wartość obie, co oznacza, że jeśli użytkownik sam nie określi wartości tego argumentu, funkcja milcząco przyjmie wartość obie. Użycie wartości domyślnej oznacza w tym przypadku, iż funkcja będzie „zachowywać się” tak jak poprzednio, nawet jeśli nie określimy żadnej wartości argumentu which.

Po drugie, wynik, jaki zwraca nasza funkcja zależy od tego, jaką wartość podstawimy pod argument which. Może on przybrać tylko jedną z trzech wartości: główna przekątna, indeks rozbieżności i obie. Mówiąc inaczej, jeśli argument which przybierze jedną z wymienionych trzech wartości, funkcja zwróci, odpowiednio, odsetek na głównej przekątnej, indeks rozbieżności lub obie te wartości. W przeciwnym przypadku otrzymamy komunikat o błędzie.

Zachowanie naszej funkcji po wprowadzeniu wszystkich tych zmian ilustrują przykłady poniżej:

mdi(tab_1)
##    główna przekątna indeks rozbieżności 
##                0.59                0.13
mdi(tab_1, which = "obie")
##    główna przekątna indeks rozbieżności 
##                0.59                0.13
mdi(tab_1, which = "główna przekątna")
## główna przekątna 
##             0.59
mdi(tab_1, which = "indeks rozbieżności")
## indeks rozbieżności 
##                0.13
mdi(tab_1[, -1], which = "indeks rozbieżności")
## Uwaga! Tylko kwadratowe macierze!
mdi(tab_1, which = "coś tam, coś tam")
## Uwaga! argument which wymaga jednej z trzech wartości: 'obie', 'indeks rozbieżności', 'główna przekątna'.
mdi(tab_1[, -1], which = "coś tam, coś tam")
## Uwaga! argument which wymaga jednej z trzech wartości: 'obie', 'indeks rozbieżności', 'główna przekątna'.

W tym miejscu się zatrzymamy, choć trzeba mieć świadomość, że niniejsza notatka stanowi w najlepszym zaledwie króciutkie wprowadzenie do wprowadzenia do tematu. Do problemu tworzenia własnych funkcji i ich zastosowania w analizach będę zapewne wracał jeszcze niejeden raz na tych stronach. Póki co, będę wdzięczny za uważną lekturę tych przykładów i ewentualne pytania.