---
title: "Wstęp do testów nieparametrycznych"
output:
  pdf_document: default
  html_document: default
---

# Wstęp do testów nieparametrycznych

### Parametryczne vs nieparametryczne

Większość poznanych przez Państwa procedur statystycznych obejmuje estymację parametrów rozkładu populacji, z której pochodziła próba lub próby. Procedury te opierają na pewnych założeniach dotyczących rozkładu populacji. Na przykład, jednym z założeń testu *t* jest to, że populacja, z której pochodzą próby (lub próba) ma rozkład normalny. Takie testy (opierające się na estymacji parametrów i/lub założeniach dotyczących rozkładu) nazywać będziemy testami **parametrycznymi**. W odróżnieniu od nich testy **nieparametryczne**, które nie opierają się na estymacji parametrów rozkładu i nie zakładają niczego o rozkładzie populacji, z której pochodzi próba lub próby.

### Kiedy użyć testu nieparametrycznego?

Przede wszystkim powinniśmy się zastanowić, czy nasze dane spełniają założenia testu parametrycznego, którego chcemy użyć. Jeżeli chcemy użyć testu $t$, to musimy zadać sobie pytanie - na ile uprawnieni jesteśmy, żeby zakładac, że rozkład danej zmiennej w populacji jest normalny. 

#### Kiedy możemy rozsądnie założyć, że rozkład w populacji nie jest normalny?
- kiedy nasza zmienna jest na skali porządkowej
- kiedy mamy jakieś obserwacje odstające (których nie możemy się pozbyć)
- kiedy widać to na jednym z wykresów diagnostycznych (histogram, wykres kwantyl-kwantyl)

#### Kiedy jeszcze możemy użyć testów nieparametryczny?
- kiedy bardziej niż średnia interesuje nas mediana
- kiedy mamy bardzo małe grupy

#### Kiedy nie jest rozsądnie założyć, że rozkład populacji nie jest normalny?
- kiedy mamy dużą próbę i jakiś test na normalność odrzuci nam hipotezę, że dane pochodza z rozkładu normalnego, ale w rzeczywistości rozbieżność nie jest duża (testy parametryczne sa dość odporne na lekką skośność itp.)

## Testy na rangach

### Test sumy rang Wilcoxona

Załóżmy, że chcielibyśmy przetestować hipotezę zgodnie z którą długofalowe bieganie redukuje ciśnienie krwi. Aby przetestować tę hipotezę zmierzyliśmy ciśnienie krwi naszym uczestnikom a nastepnie kazaliśmy im systematycznie biegać przez 6 miesięcy, po których ponownie zmierzyliśmy ciśnienie krwi.

Naszą zmienną zależną jest w tym wypadku *zmiana ciśnienia krwi po 6 miesiącach*. Jeżeli bieganie rzeczywiście obniża ciśnienie krwi, to oczekujemy że większość badanych uzyska niższy wynik przy drugim pomiarze. Oczekujemy więc dodatniej różnicy między pomiarami. Oczekujemy również że wśród tych, u których ciśnienie krwi zwiększyło się, zmiana ta nie będzie duża.

Z drugiej strony, jeżeli bieganie jest nieefektywne, to oczekujemy że około połowa różnic będzie dodatnia, a połowa ujemna, a różnice dodatnie będą mniej-więcej takie same jak te ujemne, jeśli chodzi o ich odległość od zera.

Załóżmy, że opisany powyżej eksperyment dał następujące wyniki:

```{r}
przed <- c(130, 170, 125, 170, 130, 130, 125, 160)
po <- c(120, 163, 120, 135, 143, 136, 124, 120)
roznica <- przed - po

print(przed)
print(po)
print(roznica)
```

```{r}
hist(roznica)
```

```{r}
qqnorm(roznica)
qqline(roznica)
```

Przyjmijmy, że nie możemy (albo nie chcemy!) zakładać, że nasze obserwacje pochodzą z rozkładu normalnego. Co możemy w takiej sytuacji zrobić? Czy jedyne co nam pozostaje, to poddać się? Nie! Tak jak powiedzieliśmy sobie wcześniej, jeśli bieganie nie ma wpływu na ciśnienie, to oczekujemy (odrobinę upraszczając), że połowa naszych różnic będzie ujemna, a połowa będzie dodatnia. Mówiąc bardziej precyzyjnie, nasza hipoteza zerowa głosi, że różnice będzie symetrycznie rozłożone wokół zera, hipoteza alternatywna (którą chcemy potwierdzić!) mówi zaś, że nie będą symetrycznie rozłożone wokół zera.

Na początek wykluczymy z naszej próby wszystkie te obserwacje, w których różnica wynosiła dokładnie zero. W naszej próbie nie ma takich obseracji, ale przeprowadzając test Wilcoxona należy to zrobić.

```{r}
print(roznica)
roznica <- roznica[roznica != 0]
```

Następnie uporządkujemy nasze obserwacje od tych, gdzie wartość bezwzględna różnic jest najmniejsza do tej, gdzie jest największej.

```{r}
sorted_roznica <- roznica[order(abs(roznica))]
print(sorted_roznica)
```

Następnie przypiszemy każdej obserwacji 1 jeśli różnica jest dodatnia i -1 jeśli różnica jest ujemna (jest to funkcja *signum* nazywana także w skrócie $sgn$).

```{r}
sgn_roznica <- ifelse(sorted_roznica > 0, 1, -1)
print(sgn_roznica)
```

Następnie każdej naszej obserwacji przypiszemy rangę od $1$ do $n$ gdzie $n$ to liczba obserwacji.

```{r}
ranks <- 1:length(sorted_roznica)
print(ranks)
```

Naszą statystyką testową (którą nazwiemy $V$) będzie suma iloczynów *signum* różnicy i jej rangi dla tych obserwacji, dla których signum jest dodatnie.

```{r}
V <- sum(sgn_roznica[sgn_roznica > 0] * ranks[sgn_roznica > 0])
print(V)
```

Jaki rozkład ma tak zdefiniowana zmienna losowa? Sprawa nie jest prosta - dla małych liczebności prób rozkład ten możemy wyliczyć (ale jest to dość zasobochłonne!) dla dużych ($n>50$) korzysta się czasami z przyblizenia z rozkładu normalnego. My jednak w ten problem nie będziemy się wgłębiać i posłużymy się odpowiednią funkcją wbudowaną w R w celu obliczenia wartości $p$:

```{r}
psignrank(V-1, length(roznica), lower.tail = F)
```

Nic oczywiście nie stoi na przeszkodzie, aby posłużyć się tutaj funkcją `wilcox.test` wbudowaną w R:

```{r}
wilcox.test(przed, po, paired = TRUE)
wilcox.test(przed, po, paired = TRUE, alternative = 'greater')
```

**UWAGA** W literaturze spotkają się Państwo z nieco innym sposobem obliczania statystyki testowej w tym przypadku (i inną jej nazwą, $W$). Proszę uwierzyć, że te sposoby są ekwiwalentne (oczywiście jeśli korzystamy z właściwego rozkładu!)

Niestety, okazuje się, że 6 miesięcy biegania nie sprawi, że nagle staniemy się okazami zdrowia! A szkoda...

### Test U Manna-Whitneya/Test sumy rank Wilcoxona dla prób niezależnych

Rozważmy następujące dane dotyczące liczby stresujących wydarzeń w ostatnim czasie wśród pacjentów oddziału kardiologicznego i ortopedycznego. Wiemy, że stresującym zdarzeniom w życiu (małżeństwo, nowa praca, śmierć partnera, urodzenie dziecka) często towarzyszą choroby i jest rozsądnie oczekiwać, że pacjenci oddziału kardiologicznego przeważnie doświadczyli większej liczby stresujących wydarzeń, niż pacjenci oddziału ortopedycznego. Mamy pewne podejrzenia, że wyniki na skali stresu mogą nie mieć symetrycznego rozkładu w populacji (szczególnie w przypadku pacjentów kardiologicznych, jeżeli nasza hipoteza potwierdziłaby się) dlatego użyjemy testu nieparametrycznego.


Tym razem nie będziemy przeprowadzać obliczeń ręcznie - są dość podobne do obliczeń wykonanych w poprzednim zadaniu. Zainteresowani doczytają w jaki sposób przeprowadzić ten test "na papierze". My posłużymy się funkcją `wilcox.test`.

```{r}
orto <- c(10, 12, 17, 13, 19, 20)
kardio <- c(30, 26, 25, 33, 18, 27)
data_ex <- data.frame(stres = c(orto, kardio), 
                      grupa = c(rep('orto', length(orto)), rep('kardio', length(kardio))))

wilcox.test(stres ~ grupa, alternative = "greater", data = data_ex)
boxplot(stres ~ grupa, data = data_ex, main = "orto vs kardio")
```

## Zgodność sędziów

### Kappa Cohena

Załóżmy, że poprosiliśmy sędziego ze znacznym doświadczeniem klinicznym o przeprowadzenie wywiadów z 30 nastolatkami i zakwalifikowanie każdego z nich jako przejawiających: (1) brak problemów behawioralnych, (2) internalizację problemów behawioralnych, (3) eksternalizację problemów behawioralnych.

Każdy, kto recenzowałby naszą pracę zastanawiałby się nad wiarygodnością naszych pomiarów - skąd wiemy, że sędzia jest w tym co robi lepszy, niż gdyby po prostu zgadywał na chybił-trafił? Aby wykluczyć tę możliwość, poprosiliśmy drugiego sędziego o porównywalnym wykształcenie oraz doświadczeniu zawodowym o przeprowadzenie niezależnej klasyfikacji na podstawie transkrypcji tych samych wywiadów. Po zebraniu danych, to znaczy przeprowadzeniu wywiadów oraz ich ocenie przez dwóch sędziów, utworzyliśmy tabelę krzyżową pokazującą zgodę (oraz jej brak) między naszymi oceniającymi.


```{r}
# Potrzebny nam jest pakiet `irr` (interrater reliability)
#install.packages('irr')

judge_1 = c(rep('brak problemu',15), rep('internalizacja', 2), rep('eksternalizacja', 3),
            rep('brak problemu', 1), rep('internalizacja', 3), rep('eksternalizacja', 2), 
            rep('brak problemu', 0), rep('internalizacja', 1), rep('eksternalizacja', 3))
judge_2 = c(rep('brak problemu',20), rep('internalizacja', 6), rep('eksternalizacja',4))

kategoryzacja = data.frame(przypadek = 1:30, sedzia_1 = judge_1, sedzia_2 = judge_2)
kategoryzacja
table(kategoryzacja$sedzia_1, kategoryzacja$sedzia_2)
```

Moglibyśmy po prostu policzyć, w jakiej liczbie przypadków nasi sędziowie się zgadzali i zobaczyć, czy jest to więcej niz $50\%$ przypadków. Byłoby to jednak złe podejście. Gdyby np. każdy z sędziów losowo w 80% przypadków stawiał określoną diagnozę, to wówczas mielibyśmy dość duzą szansę na to, że między nimi wystąpi dość duża zgoda, ale będzie to wciąż wynik przypadku a nie faktycznej zbieżności opinii. 

```{r}
library(irr)
agree(kategoryzacja[, c("sedzia_1", "sedzia_2")])
```

Sposobem na poradzenie sobie z tym jest obliczenie wartosci Kappa ($\kappa$) Cohena. Jest to zgodność między sędziamy z poprawką na to, że jej część moze być efektem przypadku. Możemy nawet potraktować ją jako statystykę testową i sprawdzić jakie jest jej prawdopodobieństwo przy hipotezie zerowej ($\kappa = 0$, czyli cała zaobserwowana zgodność jest wynikiem przypadku) ale proszę pamiętać, ze nawet małę wartosci będą statystycznie istotne.

W tym artykule znajdą Państwo pewne reguły interpretacji Kappy dla nauk biomedycznych (a podobne artykuły dotyczące innych dziedzin łatwo znaleźć w internecie):

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3900052/

```{r}
kappa2(kategoryzacja[, c("sedzia_1", "sedzia_2")])
```

### Kappa Fleissa

A co by się stalo, gdybyśmy mieli nie 2 a 3 sędziów? Musielibyśmy wówczas użyć innej miary - zamiast Kappy Cohena użyjemy wtedy Kappy Fleissa. Jej interpretacja jest analogiczna.

```{r}
judge_3 = c(rep('brak problemu',17), rep('internalizacja', 8), rep('eksternalizacja',5))
kategoryzacja2 = data.frame(sedzia_1 = judge_1, sedzia_2 = judge_2, sedzia_3 = judge_3)
kappam.fleiss(kategoryzacja2)
```

## Nieparametryczna "analiza wariancji"

### Test Kruskala-Wallisa

Jest to nieparametryczny odpowiednik analizy wariancji, z tą różnicą, że zamiast na wartościach pracujemy na rangach, tak jak w przypadku testu sumy rang Wilcoxona (test Kruskala-Wallisa jest jego uogólnieniem).

Załóżmy, że chcielibyśmy dowiedzieć się czy w którymś miesiacu roku poziom ozonu jest większy, niż w innych. W tym celu od maja do września przeprowadzaliśmy pomiary stężenia ozonu w powietrzu. Załóżmy również, że mamy dobre powody aby uznać, że rozkład, z którego pochodzą nasze obserwacje odbiega od rozkładu normalnego. Dodatkowo mamy pewne powody, żeby mniejszą wagę przyznawać obserwacjom skrajnym - to, że w jakimś miesiącu jeden czy dwa dni charakteryzowały się dużym stężeniem ozonu nie powinno mieć duzego znaczenia.

W tym przypadku hipoteza zerowa głosi, że obserwacje w każdym miesiącu pochodzą z tego samego rozkładu.

```{r}
boxplot(Ozone ~ Month, data = airquality)
kruskal.test(Ozone ~ Month, data = airquality)
```

Sam test Kruskalla-Wallisa nie powie nam które miesiące statystycznie odbiegają od innych. Mówi nam tylko, że przynajmniej jedna próba pochodzi z innego rozkładu niż inna próba. Jak dowiedzieć się które to próby i które miesiące są istotnie różne od których? Musimy posłużyć się testem *post-hoc*. W przypadku testu Kruskala-Wallisa musimy wybrać taki, który również operuje na rangach. W tym przypadku odpowiedni jest test Dunna.

```{r}
#install.packages('dunn.test')
library(dunn.test)
dunn.test(airquality$Ozone, airquality$Month)
```

W tabeli w każdej komórce mamy (na górze) wartosć statystyki testowej oraz (na dole) wartość *p*. Z tebeli możemy, analizując przecięcia kolumn i wierszy, odczytać, że statystycznie istotną różnicę możemy zaobserwować między 5 i 7 miesiącem, między 5 i 8, między 6 i 7, 6 i 8, 7 i 9 oraz 8 i 9. 

## Wariacje na temat $\chi^2$

### Test McNemara

Wyobraźmy sobie, że chcemy sprawdzić, czy uczestnictwo w kursie "Statystyka z R" zwiększa sympatię do statystyki. W tym celu przeprowadziliśmy badanie - spytaliśmy grupę 64 studentów kognitywistyki przed i po kursie o ich uczucia do statystyki.

Moglibyśmy skorzystać ze zwykłego testu niezależności $\chi^2$, ale nie uwzględnilibyśmy wtedy faktu, że pomiary można połączyć w pary. Zrobilibyśmy to wtedy tak:

```{r}
statystyka <- matrix(c(36, 28, 45, 19),
       nrow = 2,
       dimnames = list("Przed statystyką" = c("Przed statystyką", "Po statystyce"),
                       "Po statystyce" = c("Lubię", "Nie lubię")), byrow = TRUE)
statystyka
chisq.test(statystyka)
```

Możemy również sparować obserwacje i postawić trochę inną hipotezę - że częściej zmiana była w kierunki nie lubię -> lubię niż lubię -> nie lubię.

Hipoteza zerowa stwierdza więc, że prawdopodobieństwo zaklasyfikowania do komórki \[i,j\] i \[j,i\] jest takie samo.

```{r}
statystyka <- matrix(c(33, 12, 3, 16),
       nrow = 2,
       dimnames = list("Przed statystyką" = c("Lubię", "Nie lubię"),
                       "Po statystyce" = c("Lubię", "Nie lubię")))
statystyka

mcnemar.test(statystyka)
```

Jak widzimy, w przypadku testu niezależności $\chi^2$ różnica nie jest statystycznie istotna, a w przypadku testu MacNemara jest. Jest tak, ponieważ używajac tego drugiego, używamy dodatkowej informacji o obserwacjach. W praktyce skutkuje to większą mocą testu. Warto więc użyć testu McNemara, jeśli nasz schemat eksperymentalny zakłada dwa pomiary!

## Regresja logistyczna (od strony R)

### Statystyczna istotność współczynników w regresji logistycznej

Większość z państwa będzie uczyć się regresji logistycznej na zajęciach ze Statystyki II. Tutaj chciałbym Państwu pokazać tylko przykład, jak można za pomocą R dopasować odpowiedni model do danych i jak odczytać wyniki odpowiednich testów dotyczących współczynników.

Załóżmy, że chcielibyśmy dowiedzieć się, które czynniki miały wpływ na to, czy w katastrofie Titanica ktoś przeżył, czy nie. Tak się składa, że dysponujemy dość dokładnymi danymi z katastrofy (pakiet `titanic`) i możemy spróbować odpowiedzieć na to pytanie. Poniżej znajduje się objaśnienie poszczególnych zmiennych.

**Data Dictionary**
```
Variable	Definition	Key
survival	Survival	0 = No, 1 = Yes
pclass	Ticket class	1 = 1st, 2 = 2nd, 3 = 3rd
sex	Sex	
Age	Age in years	
sibsp	# of siblings / spouses aboard the Titanic	
parch	# of parents / children aboard the Titanic	
ticket	Ticket number	
fare	Passenger fare	
cabin	Cabin number	
embarked	Port of Embarkation	C = Cherbourg, Q = Queenstown, S = Southampton
Variable Notes
pclass: A proxy for socio-economic status (SES)
1st = Upper
2nd = Middle
3rd = Lower
```

```{r}
#install.packages('titanic')
library(titanic)
head(titanic_train)
model = glm(Survived ~ Sex + Age + Pclass + Fare + Embarked + SibSp,
            family = 'binomial',
            data = titanic_train)
summary(model)
```