Dzisiejszą notkę poświęcam testowi sumy rang Wilcoxona (Wilcoxon rank sum test). Dotychczas nie korzystałem z niego prawie wcale, ale sięgnałem po niego przy okazji analizy danych z eksperymentu, który przeprowadziłem niedawno z Kingą Wysieńską-Di Carlo. Analiza danych z eksperymentu jest w większości przypadków znacznie prostsza niż analiza danych z badania sondażowego i na ogół wystarczy skorzystać z testu t lub ANOVA, tyle tylko, że oba te testy wymagają założenia o normalności rozkładu zmiennej zależnej. W naszym eksperymencie zmienna zależna ma rozkład daleki od normalnego, jak widać na ząłączonym obrazku:

Uczestnicy naszego eksperymentu rozgrywali pewien wariant dylematu więźnia. Na początku gry badani otrzymywali od nas 10 PLN i musieli zdecydować, jaką część tej kwoty zostawić sobie, jaką zaś oddać partnerowi, o którym sądzili, że został losowo wybrany spośród innych uczestników eksperymentu. Jak wyjaśnialiśmy w instrukcjach przed rozpoczęciem gry, partner także otrzymywał od nas 10 PLN i miał podjąć taką samą decyzję. Każda złotówka oddana przez nas partnerowi była przez nas mnożona przez 2; podobnie każda złotówka otrzymana od partnera. Oddając partnerowi złotówkę, badany zostawał z 9 PLN, zaś jego partner zyskiwał 2 PLN i tym samym miał 12 PLN. Po tym „transferze” razem mieli o złotówkę więcej niż przed nim. Gdyby partner odwdzięczył się badanemu złotówką ze swojego „konta”, zostałby z 11 PLN, a badany zyskałby 2 PLN i tym samym również miałby 11 PLN. W efekcie, łącznie mieliby o 2 PLN więcej niż gdyby nie przekazali sobie ani grosza. Takie wzajemne transfery były więc dla graczy korzystne, bo powiększały ich stan posiadania; gdyby każdy przekazał drugiemu całą kwotę otrzymaną od nas, obaj zakończyliby grę z 20 PLN. Kłopot polegała jednak na tym, że wielkość wypłaty badanego zależała nie tylko do tego, co sam zrobił, lecz również od tego, jak zachował się partner. Oddając pieniądze partnerowi badany mógł dużo zyskać, gdyby partner zachował się podobnie, lub dużo stracić, gdyby partner zachował wszystkie swoje środki dla siebie. Z tego punktu widzenia, kwota oddana partnerowi stanowi miarę zaufania. Jak widzimy na rysunku powyżej, ponad 40% badanych przekazało partnerowi wszystkie swoje środki. Ponad 20% przekazało połowę. Rozkład jest mocno więc asymetryczny, dlatego użycie testu \(t\) do jego analizy jest nieuprawnione.

Istotne jest przy tym to, że badani różnili się pod względem statusu; różnice statusu między badanymi wywołano we wcześniejszych częściach eksperymentu. Mniejsza o szczegóły; dla naszych celów ważne jest jedynie to, że przed rozpoczęciem dylematu więźnia znali status swój i partnera. Przewidywaliśmy, że wysoki status będzie sprzyjał zaufaniu w tym znaczeniu, że gracz o wysokim statusie będzie uważany za bardziej „wiarygodnego”. Aby sprawdzić to przewidywanie, chcemy porównać zachowanie badanego o wysokim statusie, którego partner również ma wysoki status (w skrócie: para "High-High", albo "H-H"), z badanym, który sam ma wysoki status, zaś jego partner — niski (para "H-L"). Co więcej, uwzględniam tu tyle dane dotyczące gry sekwencyjnej, w której gracze podejmowali decyzje po kolei. Jeśli nasze przewidywanie jest słuszne, można oczekiwać, że badani w parach pierwszego rodzaju będą przekazywali swoim partnerom większe kwoty, niż badani w parach drugiego rodzaju.

W przykładach poniżej wykorzystuję następujące zmienne:

Rysunek poniżej porównuje rozkłady tej ostatniej zmiennej ze względu na status partnera.

Test sumy rang Wilcoxona

Test sumy rang Wilcoxona pozwala na porównanie dwóch rozkładów bez konieczności przyjmowania założeń na temat rodzaju tych rozkładów. W szczególnym przypadku, gdy oba rozkłady są identyczne co do kształtu i rozproszenia, a różnią się jedynie wartością mediany, test Wilcoxona jest w istocie testem różnicy median. Mówiąc dokładniej, hipoteza zerowa w tym szczególnym przypadku orzeka, iż mediany populacyjne obu rozkładów są jednakowe. Test sumy rang sprawdza tę hipotezę względem hipotezy alternatywnej, iż mediany populacyjne się różnią.

W przypadku ogólnym, gdy porównywane rozkłady różnią się kształem lub rozproszeniem, hipoteza zerowa w teście Wilcoxona ma inną postać. Niech \(X\) oznacza wartości analizowanej zmiennej w pierwszej grupie, zaś \(Y\) — wartości tej zmiennej w grupie drugiej. Porównując każdą wartość \(X\) z każdą wartością \(Y\) możemy określić, jak często zachodzi \(X>Y\). Hipoteza zerowa w teście Wilcoxona orzeka, iż \[P(X>Y)=P(Y>X)\]

Analiza

Wykonanie testu Wilcoxona w R jest proste. Zaczniemy od porangowania wartości zmiennej value w naszym zbiorze:

d <- d %>% mutate(r_value = rank(value))

Statystyka testowa, oznaczmy ją przez \(W\), jest po prostu sumą rang w mniejszej z dwóch grup. W naszym przypadku jest to grupa "L-H", jak pokazuje tabela w okienku poniżej:

test <- d %>% group_by(p_status) %>% summarise(n = n(), w = sum(r_value))
## Source: local data frame [2 x 3]
## 
##   p_status     n     w
##      (chr) (int) (dbl)
## 1        H    17   333
## 2        L    14   163

Mamy zatem:

sum(d$r_value[d$p_status == "L"])
## [1] 163

Aby ocenić, czy na podstawie tej informacji można odrzucić hipotezę zerową, możemy posłużyć się jedną z dwóch metod: „dokładną” lub opartą na przybliżeniu normalnym statystyki \(W\). Przyjrzymy się najpierw tej drugiej metodzie.

Jeśli hipoteza zerowa jest prawdziwa, statystyka \(W\) ma wartość oczekiwaną równą: \[E(W) = N_1 \bar{R}, \] gdzie \(N_1\) to liczebność mniejszej z dwóch grup, zaś \(\bar{R}\) to średnia rang w całej próbie. Ponieważ nasza analiza obejmuje 31 badanych, średnia rang w próbie wynosi: \[\bar{R} = \frac{1 + 31}{2} = 16,\] a zatem \[E(W) = N_1 \bar{R} = 224\] Odchylenie standardowe statystyki \(W\) wynosi natomiast: \[D(W) = \sigma_{R} \sqrt{\frac{N_1 N_2}{N_1 + N_2}} = 8,79\cdot\sqrt{\frac{14\cdot17}{31}} = 24,3\] gdzie \(\sigma_{R}\) oznacza odchylenie standardowe rang. Standaryzacja statystyki \(W\) zgodnie z uzyskanymi wartościami daje więc:

# Statystyka W
W <- sum(d$r_value[d$p_status == "L"])
## [1] 163
# Wartość oczekiwana W
mW <- sum(d$p_status == "L") * (1 + nrow(d))/2
## [1] 224
# Odchylenie standardowe W
sW <- sd(d$r_value) * sqrt(sum(d$p_status == "L") * sum(d$p_status == 
    "H")/nrow(d))
## [1] 24.3
# Standaryzacja W
(W - 0.5 - mW)/sW
## [1] -2.53

Wielkość \(0.5\) w formule powyżej to tzw. poprawka na ciągłość. Wartość \(p\) obliczamy następująco:

pnorm((W - 0.5 - mW)/sW, lower.tail = TRUE)
## [1] 0.00577

Innymi słowy, uzyskanie statystyki \(W\) równej 163 lub mniej jest mało prawdopodobne, jeśli hipoteza zerowa jest poprawna. Innymi słowy, mamy wystarczająco dużo dowodów, aby odrzucić hipotezę zerową o jednakowości rozkładów. Mówiąc jeszcze inaczej, możemy odrzucić hipotezę zerową, zgodnie z którą badani o wysokim statusie ufają swoim partnerom o wysokim i niskim statusie w jednakowym stopniu.

Metoda dokładna

Przybliżenie normalne powinno się jednak stosować wtedy, gdy próba jest odpowiednio duża i nie zawiera węzłów. Jeśli próba jest mała lub na wartościach zmiennej zależnej tworzy się wiele węzłów, należy skorzystać z metody dokładnej, tzn. metody, która pozwala wyznaczyć dokładną wartość \(p\) za pomocą tzw. testu permutacyjnego. Test ten polega na zastąpieniu wartości zmiennej zależnej ich losową permutacją. Mówiąc inaczej, oryginalne wartości zmiennej zastępujemy ich losową permutacją, następnie obliczamy statystykę \(W\) według podanej wyżej procedury, po czym wartości zmiennej zależnej ponownie zastępujemy kolejną permutacją losową i znowu obliczamy statystykę tetstową, itd., aż do ,,wyczerpania’’ wszystkich możliwych permutacji wartości zmiennej zależnej. Z tego też względu, opisywaną metodą daje się sensownie stosować jedynie dla niezbyt wielkich prób.

Ponieważ w teście sumy rang interesują nas (a) rangi w (b) mniejszej z dwóch grup, wykonanie testu permutacyjnego sprowadza się w istocie do wyznaczenia zbioru 14-elementowych (bo tyle jest jednostek w mniejszej z dwóch grup w naszym badaniu) kombinacji rang i obliczeniu sumy rang w każdej z uzyskanych. Posłużymy się w tym celu funkcją combn w R:

rank_perms <- combn(x = d$r_value, m = 14, FUN = sum)

Jak pamiętamy, statystyka \(W\) wyznaczona dla naszych danych jest równa 163. Wektor rank_perms zawiera wartości statystyki \(W\) dla wszystkich losowych permutacji zmiennej zależnej. Łączna liczba takich permutacji wynosi 265,182,525. Ponieważ, zgodnie z podanym wyżej przewidywaniem, oczekujemy, że zaufanie badanego do partnera w parach "H-L" będzie niższe niż w parach "H-H", pytanie, jakie sobie stawiamy brzmi: w ilu spośród tych przypadków, statystyka \(W\) przyjmuje wartość 163 lub mniejszą? 

mean(rank_perms <= sum(d$r_value[d$p_status == "L"]))
## [1] 0.00579

Zielona linia na rysunku powyżej oddziela wartości \(W\) nie większe od 163 do wartości większych niż 163. Jak widzimy, tych pierwszych jest daleko mniej. Wartość wyświetlona w okienku powyżej informuje nas, że permutacje, dla których zachodzi \(W \leq 163\) stanowią 0,00579 wszystkich wartości. Mówiąc inaczej, prawdopodobieństwo uzyskania wartości \(W\) nie większej niż 163, jeśli hipoteza zerowa jest poprawna wynosi mniej niż 0,01. Jest do dokładna wartość \(p\) w jednostronnym teście permutacyjnym.

Funkcje wbudowane

W przykładach powyżej pokazałem, jak wyznaczyć wynik testu Wilcoxona krok po kroku, bo uważam, że dobrze jest rozumieć, jak wykonuje się testy statystyczne. W praktyce jednak rzadko korzysta się z takiego „ręcznego” wyliczania statystyki \(W\) i odpowiadającej jej wartości \(p\), zwłaszcza, że test permutacyjny według pokazanego wyżej schematu jest dość czasochłonny nawet przy niewielkiach próbach. Szczęśliwie, ustnieje kilka funkcji w R, które zwracają wynik testu Wilcoxona.

Zacznijmy od funkcji wilcox.test(); znajduje się ona w pakiecie stats w bazowej dystrybucji R i jej użycie nie wymaga instalowania ani ładowania dodatkowych pakietów. Ponieważ przewidujemy, że zaufanie w parach "H-H" będzie większe niż w parach "H-L" wybieramy alternative = "greater:

wilcox.test(value ~ factor(p_status), data = d, alternative = "greater", 
    exact = FALSE)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  value by factor(p_status)
## W = 200, p-value = 0.006
## alternative hypothesis: true location shift is greater than 0

Statystyka \(W\) zwracana przez funkcję wilcox.test() różni się od wartości obliczonej powyżej, co wynika zapewne stąd, iż w literaturze podaje się różne sposoby wyznaczania tej statystyki, ale wartość \(p\) jest (w zaokrągleniu) taka sama.

Funkcja wilcox.test() podaje wartość \(p\) w oparciu o przybliżenie normalne. Wartość dokładną można uzyskać sięgając po funkcję wilcox_test() z pakietu coin. Argument distribution == "exact" wskazuje, iż chodzi nam o wartość dokładną:

library(coin)
wilcox_test(value ~ factor(p_status), data = d, distribution = "exact", 
    alternative = "greater")
## 
##  Exact Wilcoxon-Mann-Whitney Test
## 
## data:  value by factor(p_status) (H, L)
## Z = 3, p-value = 0.006
## alternative hypothesis: true mu is greater than 0
Informacja

Opisany wyżej eksperyment został przeprowadzony dzięki grantowi Narodowego Centrum Nauki pt. Dylematy kooperacji w warunkach zróżnicowania statusowego.