# Non--parametric tests {#np-tests}

So far, all of the statistical tests we've encountered rely on the same underlying logic. They try to answer the question:

**What is the probability of observing data like ours, or more extreme, if $H_0$, the null hypothesis, were true?**^[Recall that the term _null distribution_ simply means the distribution of the test statistic (i.e. $t$, $F$ or $\chi^2$ distribution with the appropriate degrees of freedom) under the premise of a true $H_0$.]

To answer this question, we calculated our observed test statistic $t$, $F$ or $\chi^2$ and placed it on the corresponding distribution (the null distribution of the test) to see how far out it is. To be able to do this, we obviously need to known the shape of the null distribution.  How do we know that?!

So far, the null distribution that we used are based on a set of assumptions which are concerned with either the parent distribution of the variable under consideration or the sample size. When the assumptions are met, it is possible to mathematically derive the null distributions (we won't do that here!). So what do we do when the assumptions are not met? The remedy to this situation comes in form of a set of statistical tests that relies on a relaxed set of assumptions -- the non-parametric tests.

## Assumptions

Assumptions are an important and ubiquitous aspect in classiscal Statistics^[The omnipresence of assumptions in Statistics is often critised, but actually science is full of assumptions. The study of a model organism for example only makes remote sense under the assumption that model organism behaves in a way that is comparable to the organism of actual interest. So the world in all its complexity is too much for most studies. We need sensible simplifications!]. The validity of your results can depend heavily on whether or not the underlying assumptions of the applied statistical method are met. So it's about time we talk about assumptions for a bit. Actually assumptions are everywhere! Even in our sweater weather example.

```{exercise, name="Finite versus infinite populations"}
Previously, we determined the temperature by observing what clothes people outside were wearing. If we had a thermometer, we could just read--off the temperature from there. BUT we'd rely on some underlying assumption. Can you tell which one?
```

Using a thermometer only makes sense if it's properly calibrated. Otherwise, the conclusions we'd make would be wrong and we may end up confidently walking out of the door with a T-shirt on, although it's cold outside. If you think about it, valid assumptions are extremely important in any sort of decision making, and statistical tests are no exception. So what are the assumptions that we have made in the tests we already talked about?

For one, _all_ the inferential methods we have covered so far require independent observations. That is, the value of one observation does not influence that of any other. This is an assumption that even non--parametric methods don't relax. However, we can get rid of some assumptions that $t$--, $F$--, and  $\chi^2$--test (and many other statistical techniques) share, namely

1. The data are drawn from a Normal distribution, _and/or_
2. The sample size is large.

Notice that the above assumptions are the same as  the requirements for the Central Limit Theorem, see section \@ref(CLT). This is of course no coincidence and relatively easy to explain for the $t$ and $F$--test: when we look at the formulae of the $t$ and the $F$--statistic, we find an $\bar{X}$ in it. So both test statistics are transformations of the sample mean. If this is normally distributed (i.e. when the above assumptions hold), then it is easy to work with mathemetically, and distributions of transformations can be derived. In our cases, the distibutions that correspond to the transformations are the $t$-- and the $F$--distribution, respectively.^[For the sake of completeness: The assumptions underlying the $\chi^2$ don't include Normally--distributed observations, which makes sense since it can be applied to categorical data, but the Central Limit Theorem still plays a role in the derivation of the the null distribution, which is why the $\chi^2$--test shouldn't be applied for small samples.]

What we learn from this is very important. Both the $t$-- and the $F$--test should only be applied if you have a large sample or if your observations are Normally distributed!! Otherwise the assumed shape of the null distribution is wrong and we cannot reliably assess whether our observed value of the test statistic is unusual._ So what do we do if these requirements are not met?! As we said already, the answer is non--parametric methods.

## Non--parametric methods -- pros and cons

Non--parametric methods come in different flavours and there are a few partially overlapping definitions which do not matter much for our purpose. The important thing here is that they are not inherently based on the CLT, so we are not slaves to the quite restrictive assumptions mentioned above. This is particularly good if your data comes from a costly and time--intensive experiment and your sample size is small.

That sounds like very good news indeed, but what's the drawback?! Non--parametric methods tend to work on a "coarse scale".^[By coarse, we mean that most non--parametric tests don't even use the actual observed values, but only look at their order or the number of times they pass some threshold.] They neglect a lot of information, which makes them _less sensitive_ to detecting interesting patterns in the data than their parametric counterparts. A non-parametric test will result in a higher $p$-value given the same input.^[A fancy way of saying that non--parametric tests tend to result in larger $p$--values than their parametric counterparts is that non--parametric tests are _conservative_.] You may find colleagues using non--parametric tests as the default method of choice, because they want to avoid false positives. Our recommendation is to aim for parametric methods, and only when their assumptions are not met should you revert to their non--parametric counterparts.

To see this trade--off, let's return to our sweater weather example. Here, parametric methods can be compared to a thermometer which is very sensitive and tells us exactly what to wear... _if_ it's calibrated correctly. If for some reason we cannot assume this, we would turn to the coarser method of looking outside to see what other people are wearing which resembles the non--parametric methods. They are less sensitive but certainly the better choice if the assumptions inherent in the parametric methods are wrong!

## Overview

Before we dive into the realm of non--parametric tests, let's provide an overview of the parametric tests that we talked about so far and their non--parametric counterparts.

Table: (\#tab:nonparatypes) A list of parametric tests and their non--parametric counterparts.

| Parametric | Non-parametric |
|---|---|
| paired $t$--test (one--sample $t$--test) | sign--test, Wilcoxon signed rank test |
| two--sample $t$--test                    | Mann--Whitney U test (aka Wilcoxon rank sum test) |
| one--way ANOVA                           | Kruskall Wallis test |

| Large sample size | Small sample size |
|---|---|
| $\chi^2$--test    | Fisher's exact test |


In some books, the $\chi^2$--test is described as being non-parametric, despite requiring a largish sample size (see \@ref(largechi) which is why we treat the $\chi^2$--test/Fisher's exact test pair a little different in the table.

## The sign test

The sign test is probably the easiest non--parametric test and can be almost fully explained with material we covered in the chapter on probability. It is the non--parametric equivalent of the one--sample $t$--test and -- by extension -- the paired $t$--test. As the paired--data context is more interesting for practical purposes, let's investigate how the sign test works in this setup.

When we talked about the paired $t$--test, we considered two measurements that were made on the same Martians, namely the time to run 100 m as measured on Mars and then on Earth. To investigate whether we could reject the $H_0$ that the planet does _not_ influence the running speed of the Martians, we looked at the sample mean of the slopes in Figure _MartianPairedlines_. Based on the corresponding $t$--statistic, we concluded that it is further away from zero than we'd expect if the unknown mean population slope (i.e. the mean difference in running speed of _all_ Martians) were zero. 

Now we know that we really need to check the underlying assumptions of the $t$--test before we start using it: our sample size is not very big so the CLT won't do the trick. Also if we're honest we don't really have enough data to determine whether the slopes themselves are Normally distributed. So let's make sure we don't walk into the uncalibrated thermometer trap and use the sign test instead!

Before we start talking about the test statistic and the null distribution, the one thing that we need to be clear about is the exact form of the $H_0$. As we said, for the paired $t$--test it is that the unknown mean population difference (i.e. population mean slope) is equal to zero. The $H_0$ of the sign test is similar but uses the unknown population _median_ of the differences!^[Recall that the median is the robust alternative to the mean (cf. Section _location_) so this should make sense to you!]

There is no universally agreed on symbol for the population median (unlike $\mu$ which is used by almost everyone as the symbol for the population mean). Let's call it $M$. With this the $H_0$ and $H_1$^[Here, we are using a two--sided sign test, but the one--sided sign test also exists.] here could then be written as

$$ H_0:  M_{\text{Earth - Mars}} = 0 \quad \quad \text{vs.}  \quad \quad H_1:   M_{\text{Earth - Mars}} \neq 0. $$

Now that we know the exact test problem we are looking at, we can start thinking about a suitable test statistic and the corresponding distribution this statistic would have under $H_0$. In the case of the sign test, this is really straight forward -- all we need is to remember the definition of the median!

The median divides the dataset in halves. If $H_0$ were true, this would mean that half of the differences (Earth-Mars) are smaller and the other half are larger than zero. So the probability that a single difference is larger than zero is 0.5. This situation should ring quite a few bells, because this is exactly what we have in the context of a Binomial distribution. Let's define a positive difference as 'success'. With this we can use the number of times this happens as our test statistic $T$, we have everything we need! Under a true $H_0$ the test statistic follows a Binomial distribution with the parameters $p=0.5$ and $n=$ our sample size. 


```{r FirstNP, echo = FALSE, fig.cap="Observed differences in running time between Mars and Earth, the blue ones are negative and the red ones are positive.", fig.width = large_width*1.5, echo = FALSE}

martian$diffCol <- 0
martian$diffCol[martian$Ratio > 1] = 1
martian$diffCol <- factor(martian$diffCol)


paired3 <- ggplot(martian[martian$Site == "Site I",], aes (x = Group, y = Ratio, col = diffCol)) +
  coord_fixed(0.25, xlim = c(0.5,1.5)) +
  geom_point(alpha = 0.6, size = 4, position = position_jitter(0.3, seed = 136)) +
  scale_colour_manual(values = c(DarkBlue, DarkRed), guide = FALSE) +
  geom_hline(yintercept = 0)+
  scale_y_continuous("Observed Difference in Means", limits = c(-8, 4), breaks = seq(-8, 4, 2), expand = c(0,0)) +
    theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.y = element_line(colour = "black"),
        axis.line.y = element_line(colour = "black"),
        axis.text.y = element_text(colour = "black"),
          axis.line.x = element_line(color = "white"),
        axis.ticks.x = element_line(color = "white"),
        axis.title.x = element_text(color = "white"),
        axis.text.x = element_text(color = "white")
  )


# martian %>% 
#   select(c(Site, Speed.Mars, Speed.Earth)) %>% 
#   gather(key, value, -Site) %>% 
#   filter(Site == "Site I") %>% 
#   mutate(key = factor(key, labels = c("Mars", "Earth")),
#          obs = rep(1:(length(value)/2), 2),
#          jitt = jitter(c(rep(1, (length(value)/2)),rep(2, (length(value)/2))), .2)) -> martian.m
# 
# # Paired:
# martian.m %>% 
#   select(-Site, -jitt) %>%
#   spread(key, value) %>% 
#   mutate(diff = Mars - Earth,
#          group = 1) -> martian.m

cowplot::plot_grid(paired1, paired2, paired3, ncol = 3, align = "v", axis = "b")

```


```{r FirstNP1, fig.cap="The Binomial distribution with $n=10$ and $p=0.5$ is the null distribution of the sign test in our example. The black line shows the observed test statistic $T=2$ and the sum of the read bars corresponds to the $p$--value.", fig.width = large_width*1.5, echo = FALSE}

prob <- 0.5
size <- 10
obs <- 2

val <- dbinom(0:size, size, prob)
df.raw <- data.frame(variable = seq(0, size, 1), values = val)
df.martian <- as.data.frame(df.raw)
names(df.martian) <- c("variable", "values")

point.prob <- df.martian$values[df.martian$variable == obs]

lower.bound <- max(df.martian$variable[df.martian$values <= point.prob & df.martian$variable < obs])

df.martian$col <- NULL
df.martian$col[df.martian$values <= point.prob] <- MidRed # Red
df.martian$col[df.martian$values > point.prob] <- MidYellow # Yellow

# Total area:
sign_p_value <- sum(df.martian$values[df.martian$values <= point.prob])
# Same as:
# binom.test(2, 10)$p.value


# Draw basic yellow probability plot:
ggplot(df.martian, aes(x = variable, y = values, fill = col)) +
  geom_bar(stat = "identity") +
  # geom_vline(xintercept = obs) +
  # geom_hline(yintercept = df.martian$values[df.martian$variable == obs]) +
  scale_fill_identity(guide = FALSE) +

  scale_x_continuous(paste("Binomial Probability Distribution of\n", size," trials with a ", prob, " probability of success", sep = ""),
                     #labels = categories,
                     limit = c(-0.5, 10.5),
                     labels = seq(0,10,1),
                     breaks = seq(0,10,1),
                     expand = c(0,0)) +
  scale_y_continuous("Probability",
                     limit = c(0, 0.25),
                     labels = seq(0, 0.25, 0.05),
                     breaks = seq(0, 0.25, 0.05),
                     expand = c(0,0))  +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        axis.line.x = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )


```



For our Martian example, $n=10$, so the null distribution of $T$ is the distribution shown in fig. \@ref(fig:FirstNP2). The observed value of the test statistic here is $T=2$ (two observed differences are positive, see fig \@ref(fig:FirstNP)), so we can read--off the $p$--value in a familiar manner. In this case the p-value is `r sign_p_value`. 

This is pretty simple -- moreso than the $t$--test, so why don't we use this test all the time?! The sign test is a very good showcase for the strengths and weaknesses of non--parametric tests. On the one hand, it is strikingly simple and elegant in its logic, _but_ it doesn't make use of all the information in the data. A continuous variable (difference in speed) was "downgraded" to a binary variable (bigger than 0, yes or no). This is a clearly more coarse way to look at the data, and the associated loss of information is particularly bad for small samples. 

To understand this, remind yourselves about the Binomial distribution with $p=0.5$ and $n=3$. Here, the $p$--value of the two--sided sign test can either be 0.5 (if you observe 0 or 2 successes) or 1 (if you observe 1 success). So the $p$--value comes with very coarse steps.^[We know this seems contradictory, as we earlier said that non--parametric methods are the way to go with small samples. As we see here, even they cannot turn a small $n$ into a good thing. If you have to make an important decision, you need lots of evidence and a large sample size. Sorry... but don't shoot the messenger!]  In particular, with a 5% significance level, you won't be able to ever reject $H_0$, no matter what the actually observed values look like! 


```{exercise, name="Finite versus infinite populations"}
What does the above problem mean if you are interested in a significance level of $\alpha = 0.05$? Does this make sense?
```

```{exercise, name="Smallest sample sizes"}
What's the minimum $n$ you need to be able to potentially reject $H_0$ at the 5% significance level? Use the webapp!
```

Thankfully, the issue of coarse steps in the $p$--values that can potentially be obtained in a sign test gets better quickly with increasing $n$! And, to leave the topic of the sign test on a positive note, it can be applied to discrete or ordinal data which is something the $t$--test cannot offer!

## Wilcoxon signed--rank test

Another non--parametric equivalent for one--sample (or paired) $t$--tests is the Wilcoxon signed--rank test. In a way, this test offers a nice compromise between the simple but coarse sign test and the sensitive but full-of-assumptions $t$--test. That's because the observed difference values are not converted into binary observations as with the sign test, but rather their magnitude matters. They enter the test statistic in the form of _signed ranks_.^[Recall that we are still dealing with the paired differences between running time for 100 m on Earth and Mars, hense the term _difference values_. You can use any kind of measured values here, but this works nicely as an example since negative and positive values really do have different meaning.]

To obtain these signed ranks, the observed values are ordered according to their magnitude first -- starting with the smallest, _disregarding sign_, i.e. whether they are positive or negative doesn't matter! Let's use a tiny example with $n=3$ as a showcase in which we observed the following difference values $d_1 = 1.1$, $d_2 = -2.3$, $d_3 = 1.5$. The table below shows the order of their absolute values and the corresponding ranks.

<!-- % which can be viewed as the place of the absolute differences in the competition of being the smallest. -->

|    |   |   |   |
|--------------------------|-----|-----|-----|
| ordered absolute values  | 1.1 | 1.5 | 2.3 |
| ranks of absolute values | 1   | 2   | 3   |
| sign                     | +   | +   | -   | 

As the next step, we need to sum up the ranks with a positive sign. Let's call this sum $W_{+}$. Similarly, the sum of the negative ranks are denoted with $W_{-}$. In our example, $ W_{+} = 3$ and $W_{-} = 3$, and they contribute equally to the total sum of the ranks $W_{\text{ tot} = 6}$. Even without knowing what comes next, you can probably imagine that this seems supportive of the $H_0$, that the population median of the differences between the treatments is equal to zero and and the observed differences are not systematically below or greater than zero.^[Again, a one--sided alternative exist for the Wilcoxon signed--rank test as well.]

The test problem (i.e. $H_0$ and $H_1$) coincides with the sign test. Now, we still need a test statistic and a null distribution. As a test statistic you can either choose $W_{+}$ or $W_{-}$, here we choose $W_{+}$^[In some textbooks you may find other versions of the test statistic, like the minimum of $W_{+}$ and $W_{-}$. The test procedure, however, is always the same, so if you understand what we are saying here, you'll be set!]. What's its distribution under $H_0$?! For our tiny example, we can even derive this distribution by hand!

With $n=3$, we have three ranks available, 1, 2, and 3, and they can either all correspond to negative differences, or all to positive differences, or anything in between. There are $2^3 = 8$ possibilities combinations of 2 signs (+ and -) and 3 ranks as shown in Table \@ref(tab:WilcoxNP).


WilcoxNP All possible combinations of signed ranks and their associated point probabilities for a Wilcoxon signed--rank test when $n=3$.

| Combination | 1 | 2 | 3 | $W_{+}$ | $W_{-}$ | point probability | two-sided p-value |
|-------------|---|---|---|---------|---------|-------------------|-------------------|
| 1           | + | + | + | 6       | 0       | 0.125             | 0.25              |
| 2           | - | + | + | 5       | 1       | 0.125             | 0.50              |
| 3           | + | - | + | 4       | 2       | 0.125             | 0.75              |
| 4           | + | + | - | 3       | 3       | 0.125             | 1.00              |
| 5           | - | - | + | 3       | 3       | 0.125             | 1.00              |
| 6           | - | + | - | 2       | 4       | 0.125             | 0.75              |
| 7           | + | - | - | 1       | 5       | 0.125             | 0.50              |
| 8           | - | - | - | 0       | 6       | 0.125             | 0.25              |

Under the null hypothesis, all eight combinations are equally likely, so each of them can occur with probability 1/8. For our test statistic $W_{+}$, this means that _all_ values 0 - 6 have probability 1/8, except 3 which has a probability of 1/4 (1/8 + 1/8). This is all there is -- that's the null distribution of the test statistic. Figure \@ref(fig:WilcoxNP) shows the corresponding plot of the null distribution.

```{r WilcoxNP, fig.cap="The null distribution of the Wilcoxon signed--rank test when $n=3$.", eval = FALSE, echo = FALSE}

x <- as.character(0:6)
p <- rep(1/8,7)
p[4] <- 1/4

df <- data.frame(x,p)

ggplot(df, aes(x = x, y = p)) +
      geom_bar(stat = "identity", fill = MidYellow) +
      scale_x_discrete("$W\\_+$") +
      scale_y_continuous("probability", limits = c(0, 0.25), breaks = c(0, 0.125, 0.25)) +
      theme(panel.background = element_blank(),
            plot.background = element_blank(),
            legend.background = element_blank(),
            legend.key = element_blank(),
            axis.text = element_text(colour = "black"),
            axis.ticks.y = element_line(colour = "black"),
            axis.ticks.x = element_line(colour = "black"),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank()
      )

    # So let's look at our value:
    # wilcox.test(c(1.1, -2.3, 1.5), exact = TRUE, correct = FALSE) # V = 3, p-value = 1
    #
    #
    #
    # wilcox.test(c(1.1, 1.5, 2.3), exact = TRUE, correct = FALSE) # V = 6, p-value = 0.25
    # wilcox.test(c(-1.1, 1.5, 2.3), exact = TRUE, correct = FALSE) # V = 5, p-value = 0.5
    # wilcox.test(c(1.1, -1.5, 2.3), exact = TRUE, correct = FALSE) # V = 4, p-value = 0.75
    # wilcox.test(c(1.1, 1.5, -2.3), exact = TRUE, correct = FALSE) # V = 3, p-value = 1
    # wilcox.test(c(-1.1, -1.5, 2.3), exact = TRUE, correct = FALSE) # V = 3, p-value = 1
    # wilcox.test(c(-1.1, 1.5, -2.3), exact = TRUE, correct = FALSE) # V = 2, p-value = 0.75
    # wilcox.test(c(1.1, -2.3, -1.5), exact = TRUE, correct = FALSE) # V = 1, p-value = 0.5
    # wilcox.test(c(-1.1, -1.5, -2.3), exact = TRUE, correct = FALSE) # V = 0, p-value = 0.25

wilcox_p_value <- wilcox.test(martian$Ratio[martian$Site == "Site I"], exact = TRUE, correct = FALSE)$p.value

```

```{exercise, name="Finite versus infinite populations"}
Using the null distribution in Figure _WilcoxNP_, can you calculate the $p$-value for our example? What's the minimal $p$-value you can obtain when you use Wilcoxon signed--rank test with $n=3$?
```

<!-- % From this you see that even if all difference values are negative or all are positive, you wouldn't be able to reject the null hypothesis at $\alpha = 0.05$, but at least the situation is a lot better than the one we faced in the sign test where the minimal $p$--value that could possibly obtained with $n=3$ was 0.5. Here it is 0.25 -- and again, the situation gets better quickly with increasing $n$. -->

How does the Wilcoxon signed--rank test compare to the sign test? The sign test wouldn't make a distinction between combinations 2-4 in Table \@ref(tab:WilcoxNP), as they all have two positive and one negative sign. The same goes for combinations 5-7. The Wilcoxon signed--rank test, on the other hand, recognises that it does make a difference if only the smallest observed difference or only the largest observed difference is negative. The former points a little bit towards a positive (population) median whereas the latter points towards a (population) negative median. Therefore, the magnitude of the observed difference plays a role in the form of ranks which is why we said that the Wilcoxon signed--rank test works on a finer scale than the sign test.
    
If we apply the Wilcoxon signed--rank test to our Martian speed dataset, we get a $p$--value of _wilcoxpvalue_ Note that this value is between the $p$--values we got from the paired $t$--test _pairedTpvalue_ and the signed test _signpvalue_. This reflects how much information the tests take into account. The sign test works on a very coarse scale, so we shouldn't be surprised that the $p$--value is so high. But again, we remind you that the applied test should not be chosen according to which one gives the smallest $p$--value. If you want reproducible results, use the test whose assumptions are in line with your data situation!
    
There is one final thing worth noting here that will make a reappearance in the other non--parametric tests that we'll cover as well: the above derivation of the null distribution requires enumerating all different combinations of ranks. If $n$ is large, this could take a while (even for computers). So instead of using the _exact_ null distribution you could also use a Normal approximation of the null distribution that is simply based on the Central Limit Theorem. In this case, if the sample size is large enough so that the CLT applies, $W_{+}$ follows a Normal distribution with mean $\mu = n(n+1)/4$ and $\sigma^2 = n (n+1)(n+2)/24$.


<!-- % To demonstrate the above points, consider the observed differences $d_1 = 1.1$, $d_2 = -2.3$, $d_3 = 1.5$, $d_4 = 1.5$, $d_5 = 1.5$ ,  $d_6 = -1.1$, $d_7 = 0$ which would be dealt with as below -->
<!-- % \begin{center} -->
<!-- % \begin{tabular}{r|ccccccc} -->
<!-- % ordered absolute values & 0 & 1.1  & 1.1 & 1.5 & 1.5& 1.5 & 2.3 \\ -->
<!-- % ranks of absolute values & NA & 1.5 & 1.5 & 4 & 4&4 & 6 \\ -->
<!-- % sign & NA & + & - & + & + & + & - -->
<!-- % \end{tabular} -->
<!-- % \end{center} -->
<!-- % The average ranks of 1.5 come from (1+2)/3 and the average ranks of 4 come from (3+4+5)/3. -->
<!-- % %The sum of all ranks is 9 which is $n^* (n^*+1)/2$, where $n^* = 6$. -->
<!-- % -->

## Mann--Whitney U test (aka Wilcoxon rank sum test)

```{r echo = FALSE}

g1 <- c(-1.0, 1.5, 9.6)
g2 <- c(-1.3, -1.5, 2.0, -0.1)

wilcox_play <- wilcox.test(g1, g2)
wilcox_play <- wilcox_play$p.value

height_wilcox <- suppressWarnings(wilcox.test(Height ~ Site, martian))
height_wilcox <- height_wilcox$p.value
```

The Mann--Whitney U test is the non--parametric equivalent of the two--sample $t$--test.^[The Mann--Whitney U test was developed by two gentlemen, Mann and Whitney. Mr Wilcoxon, who we already encountered above, developed an equivalent procedure independently - which is why the Mann-Whitney U test is sometimes also called the Wilcoxon rank sum test. We call it the Mann--Whitney U test to avoid confusion with the Wilcoxon signed-rank test, above.]

It compares the unknown location parameters of two independent samples. The required input for the test are observations that come from two different groups $G_X$ and $G_Y$ with potentially different sample sizes. There are many ways to write down the considered test problem.^[For more information about the different ways to write down $H_0$ and $H_1$, see also: Morten W. Fagerland and Sandvik (2009): The Wilcoxon Mann Whitney test under scrutiny, Statistics in Medicine 28:1487--1497.] Here, we use a version that --  written as a formula -- looks a bit funny
<!-- % Here, the test problem written as a formula looks a bit funny -->
<!-- % $$ H_0 : F_X(x) = F_Y(x) \quad \quad \mbox{vs} \quad \quad H_1 :F_X(x) = F_Y(x+M).$$ -->
<!-- % In words, the $H_0$ is that the variable we measure has the exact same distribution in both groups, implying the same underlying location parameter. The $H_1$ states that the unknown population distribution in group $G_Y$ is the same as in groups $G_X$ \emph{except for a shift in location of magnitude $M$.}\marginnote[-1cm]{As an example of a location shift, you can think of the Normal distribution app where you fix a $\sigma^2$ and only change the population mean $\mu$ to $\mu + M$ In the current context, however, we are not restricted to the Normal distribution.} -->
<!-- %  -->
<!-- % If you think about it, the test problem as stated above has a strong inherent assumption which is that the underlying distribution functions differ \emph{at most} by some location shift and are exactly the same otherwise. In practice, it will be hard to verify this assumption in which you can still use the Mann--Whitney U test, but the test problem that you are investigating also looks funny but different -->

<!-- % \begin{equation} -->

$$ H_0 : P(X < Y) =  P(X \geq Y)\quad \quad \text{vs} \quad \quad H_1 :P(X < Y) \neq  P(X \geq Y).$$

<!-- % \end{equation} -->

In words, it's easy to understand though: the null hypothesis is that values from group $G_X$ have the same probability to be greater than values from $G_Y$  as they have to be smaller. The alternative hypothesis is the opposite of this, i.e. the values from one group are systematically smaller than from the other. 
<!-- % Note that this is more general than the test problem above and does not explicitly test whether the two groups differ in location. If the two underlying distributions have the same median but a different skew, the null hypothesis would also be rejected.\marginnote{This shows again that the statements that you can make about your data get coarser the fewer assumptions you can make.} -->

The actual test procedure will be quite straightforward to explain as it is very similar to the Wilcoxon rank sum test in the previous section. To see this, let's use a little toy example again where $x_1 = -1.0$, $x_2 = 1.5$, $x_3 = 9.6$ and $y_1 = -1.3$, $y_2 = -1.5$, $y_3 = 2.0$ $y_4 = -0.1$ are two samples taken from two groups $G_X$ and $G_Y$ with $n_X = 3$ and $n_Y = 4$ observations, respectively.  The first step is to rank the entire set of $n_X + n_Y$ observations irrespective of their group membership, and then to calculate the sums of the ranks from the two groups separately as show below

<!-- ordered values & -1.5 & -1.3  & -1.0 & -0.1 & 1.5& 2.0 & 9.6 \\ -->
<!-- ranks  & 1 &2 & 3 & 4 &5&6 &7 \\ -->
<!-- group & Y & Y & X & Y & X &Y &X -->

Let's denote the sums of the ranks of the two groups by $W_X$ and $W_Y$ so that $W_X = 15$ ($3 + 5 + 7$) and $W_Y=18$ (1 + 2 + 4 + 6) in our example^[Fun fact: in our example, the total sum of all ranks is $(n_X + n_Y)(n_X + n_Y+1)/2 $. This is no coincidence but a general result. Knowing this, one can calculate $W_X$ from $W_Y$ and vice versa by subtracting the other one from $(n_X + n_Y)(n_X + n_Y+1)/2 $.].
What do we expect to find if the null hypothesis were true? If neither group produces systematically larger or smaller values than the other, we expect to find a similar number of high and low ranks in each group. So if we add them up, then the group--specific sums should be about the same. Or should they?

Remeber that the two groups are allowed to have different sample sizes! So one group with many small ranks can have a similar $W$--value as a group with a few large ranks which wouldn't be very supportive of the null hypothesis. So instead of using $W_X$ and $W_Y$ directly, we need to adjust them using the group--specific sample sizes:

$$W^*_X = W_X - \frac{n_X (n_x + 1)}{2} \quad \quad \text{and} \quad \quad W^*_Y = W_Y - \frac{n_Y (n_Y + 1)}{2}$$

which for our example results in $W^*_X = 15-6=9$ and $W^*_Y=18-10=8$.
Either of them could be used as the test statistic for the Mann--Whitney U test because one can be calculated from the other. The convention is to look at the smaller value of the two which can be denoted as $W^*_{\text{min}} = \min \{W^*_X,W^*_Y\}$. 

To test our null hypothesis, we need now to know if $W^*_{\text{min}}$ we observe is unexpected under $H_0$ by comparing it to the correct null distribution. We won't go into detail here, but like in the Wilcoxon signed ranks test, the null distribution can either be derived exactly by enumerating all the potential values of our test statistic and counting how many times the unique values occur under $H_0$ or (if the sample sizes are large) using a Normal approximation of the null distribution.^[It's worth finding out, which null distribution your particular software uses!] The $p$--value in our play dataset with regard to the exact null distribution is `r wilcox_play`.

**Exercise** What's the value of  $W^*_{\text{min}}$ for the following dataset? $x_1 = -1$, $x_2 = -1.5$, $x_3 = -9.6$ and $y_1 = 1.3$, $y_2 = 1.5$, $y_3 = 2.0$ $y_4 = 0.1$

**Exercise** What's the value of  $W^*_{\text{min}}$ for the following dataset? $x_1 = -100$, $x_2 = -150$, $x_3 = -960$ and $y_1 = 130$, $y_2 = 150$, $y_3 = 200$ $y_4 = 10$

**Exercise** Compare the two results above, do they make sense? Would the $p$--values be different?

**Ties and zeros**

The procedure of the Wilcoxon signed ranks test probably raised two burning questions which we should answer before we move on to the next tests that are also based on ranks:
1. _What happens if several absolute differences have the same value?_

This problem is called ties, but luckily it isn't a real problem. The rank of all the absolute values that are the same is simply the average of the ranks they span.

2. _What happens if we observe differences of zero?_

These don't play a role, they are omitted before the whole procedure starts.
Important: This means that you have to adjust the $n$ you are working with, i.e. if 3 of your observed differences are zero, then effectively you are working with a new $n^*$ of $n-3$.


## Kruskal Wallis test

```{r echo = FALSE}

id<-(1:15)
libido<-c(3,2,1,1,4,5,2,4,2,3,7,4,5,3,6)
group <- rep(1:3,each = 5)
unique <- rep(1:5, 3)
dose<-gl(3,5, labels = c("Placebo", "Low Dose", "High Dose"))
Group<-gl(3,5, labels = c(" 1", " 2", " 3"))

viagraData<-data.frame(dose, libido, group, unique)
viagraData <- plyr::ddply(viagraData, c("dose"), transform,
                        mean = mean(libido),
                        sd = sd(libido))

KW.df <- viagraData
KW.df$rank <- rank(KW.df$libido)

rValues <- tapply(KW.df$rank, KW.df$group, mean)
# mean(KW.df$rank)
nValuesKW <- tapply(KW.df$rank, KW.df$group, length)

# sum((KW.df$rank - mean(KW.df$rank))^2)

# nominator:
nominator <- sum(nValuesKW * (rValues - mean(KW.df$rank))^2)

denominator <- sum((KW.df$rank - mean(KW.df$rank))^2)

# 14 * nominator/denominator

col.y <- "#F5B319"
col.r <- "#C42126"

#LibidoDose.lm <- lm(libido ~ dose)
LibidoDose.lm <- lm(libido ~ Group)

LibidoDose.lm.summary <- summary(LibidoDose.lm)
LibidoDose.lm.anova <- anova(LibidoDose.lm)
ANOVA.p <- LibidoDose.lm.anova$`Pr(>F)`[1]

myKW <- kruskal.test(libido ~ group, viagraData)

# Irina's example:
values.I <- c(0.1, 0.7, 0.9, -2.1,
              -0.7, 1.0, -1.2, -0.1, 3,
              1.3, 1.7, -10, 2.0, 1.3, -1.5)

groups.I <- c(1,1,1,1,
              2,2,2,2,2,
              3,3,3,3,3,3)
# rank(values.I)
ranks.I <- c(7, 8, 9, 2,
             5, 10, 4, 6, 15,
             11.5, 13, 1, 14, 11.5, 3)
# mean(ranks.I)

# sum((ranks.I - mean(ranks.I))^2)

Irina.KW <- data.frame(values.I, groups.I)

# kruskal.test(values.I ~ groups.I, Irina.KW)

```

Just like one--way ANOVA extends the two-sample $t$--test to the comparison of several groups, the Kruskal Wallis test is an extension of the Mann--Whitney U test.
<!-- % The underlying assumptions are therefore similar to those we talked about in the previous section: you don't need to assume Normally distributed data. However, if you want to investigate whether all groups have the same unknown population median, you need to assume that the shape of the population distributions in all groups is the same apart from maybe a location shift. Like in the context of the Mann--Whitney U test, the test problem needs to be modified (softened if you wish) if you cannot make this assumption. In this case, t -->
It investigates whether there is at least one group that systematically produces higher values than one of the other groups.

So how does it work?  In a nutshell, the procedure is to rank the data irrespective of the groups membership just like in the Mann--Whitney U test. Then the similarity with one--way ANOVA kicks in and the group--specific averages of the ranks as well as the total average of all ranks are calculated^[Recall that we also calculated the overall and the individual group means in ANOVA.]. Similar to the $F$--statistic, the test statistic of the Kruskal Wallis test compares the variability of the data _between_ the group to the total variability. This time, however,  we are looking at something that resembles $R^2$ (the ratio between model sum of squares and total sum of squares) and everything is done in terms of the ranks. The interpretation remains the same though: if most of the total variability is present _between_ the groups rather than _within_, then there is evidence in the data that the group--effect cannot be neglected.

Table _KWTABs_ shows the same example dataset (3 groups and 15 observations) that we used when learning about ANOVA. To understand how the ranks were calculated, look at Figure _BrokRanks_ and the margin note on page _tiesZero_.


<!-- \multicolumn{2}{c|}{Group 1}  & \multicolumn{2}{c|}{Group 2} & \multicolumn{2}{c}{Group 3} \\ -->
<!-- value & rank & value & rank & value & rank \\ -->
<!-- \hline -->
<!-- 3.0 & 7.0  & 5.0 & 12.5 & 7.0 & 15.0 \\ -->
<!-- 2.0 & 4.0  & 2.0 & 4.0  & 4.0 & 10.0 \\ -->
<!-- 1.0 & 1.5  & 4.0 & 10.0 & 5.0 & 12.5 \\ -->
<!-- 1.0 & 1.5  & 2.0 & 4.0  & 3.0 & 7.0 \\ -->
<!-- 1.0 & 10.0 & 3.0 & 7.0  & 6.0 & 14.0 \\ -->
_KWTABs_
An example dataset and associated ranks for each value. The ranks are the input in the Kruskal Wallis test.


```{r echo = FALSE}

df.ranks <- data.frame(libido,
                       group,
                       rank = rank(viagraData$libido),
                       posn = c(6,3,1,2,11,13,4,10,5,7,15,9,12,8,14),
                       fakeY = 1)
df.ranks$group <- factor(df.ranks$group)

segments.df <- data.frame(begin = c(1, 3, 6, 9, 12),
                          end =   c(2, 5, 8, 11, 13),
                          libido = 1,
                          group = 1,
                          rank = 1,
                          posn = 1,
                          fakeY = 1)
segments.df$group <- factor(segments.df$group)

plot1 <- ggplot(df.ranks, aes(x = posn, y = fakeY, col = group)) +
  geom_segment(data = segments.df, aes(x = begin, xend = end, y = posn, yend = fakeY), col = "black", size = 3) +
  geom_point(size = 6) +
  scale_x_continuous("Ordered values, without breaking ties",
                     breaks = 1:15,
                     limits = c(0,15)) +
    scale_colour_brewer(palette="Dark2") +
    theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.line.x = element_line(colour = "black"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "none"
)


plot1


plot2 <- ggplot(df.ranks, aes(x = rank, fill = group)) +
  geom_dotplot(stackgroups = TRUE, size = 2, col = NA) +
  scale_x_continuous("Ranked values, breaking ties",
                     breaks = sort(unique(df.ranks$rank)),
                     labels = sort(unique(df.ranks$rank)),
                     limits = c(0,15)) +
  scale_fill_brewer("Group", palette="Dark2") +
  theme(panel.background = element_blank(),
        legend.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.key = element_blank(),
        axis.ticks.x = element_line(colour = "black"),
        axis.line.x = element_line(colour = "black"),
        axis.text.x = element_text(colour = "black"),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "none"
)

plot2

# grid.arrange(plot1, plot2,nrow=2)


```
_BrokRanks_
The ranks for all 15 values in our dataset. _top_: Each value receives a rank according to its order among _all_ the data. Dark black lines indicate ties. _bottom_: Ties are broken by taking the average rank values.

```{r echo = FALSE}

# First just draw the individual points:
m <- ggplot(viagraData, aes(x = group, y = libido, group = factor(unique), col = factor(group))) +
  guides(colour = FALSE) +
  geom_point(position =  position_dodge(.3), size = 4) +
  scale_color_brewer(palette="Dark2") +
  scale_y_continuous("", limits = c(0,8), expand = c(0,0)) +
  scale_x_continuous("", limits = c(0.7,3.2), labels = c("Group 1", "Group 2", "Group 3"), breaks = 1:3) +
      theme(panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_line(colour = "black"),
        axis.line = element_line(colour = "black"),
#        axis.line.x = element_blank(),
        axis.text = element_text(colour = "black")
    )

m




```

_marginfigure_
The original, un-ranked data is found in three groups, each containing 5 observations.

The total average of all the ranks is $\bar{r} = `r mean(KW.df$rank)`$ and the group specific averages are $\bar{r}_1 = `r rValues[1]`$, $\bar{r}_2 = `r rValues[2]`$ and $\bar{r}_3 = `r rValues[3]`$. As stated above, the test statistic now  resembles $R^2$ from linear regression and in its general form looks like this

$$ K = (n-1) \frac{\sum_{i = 1}^g n_i (\bar{r}_i - \bar{r})^2}{\sum_{i=1}^g \sum_{j = 1}^{n_i} (r_{ij} - \bar{r})^2},$$

where $g$ denotes the number of groups and $n = \sum_{i=1}^g n_i$ is the total sample size. In our example $n = 5 + 5 + 5 = 15$. Moreover, the numerator is

$$`r nValuesKW[1]` \cdot (`r rValues[1]` - `r mean(KW.df$rank)`)^2 + `r nValuesKW[2]` \cdot (`r rValues[2]` - `r mean(KW.df$rank)`)^2 + `r nValuesKW[3]` \cdot (`r rValues[3]` - `r mean(KW.df$rank)`)^2 = `r nominator`$$

and can be interpreted as the between group variability of the ranks. If this is large, this points towards the fact that the groups are quite far apart. To judge whether they are apart enough to justify a significant result, this value need to be compared to the total variability of the ranks in the denominator, i.e. for our example

$$ (`r sort(KW.df$rank)[1]` - `r mean(KW.df$rank)`)^2 +(`r sort(KW.df$rank)[2]`-`r mean(KW.df$rank)`)^2 + \ldots +(`r sort(KW.df$rank)[14]`-`r mean(KW.df$rank)`)^2 +(`r sort(KW.df$rank)[15]`-`r mean(KW.df$rank)`)^2 = `r denominator`$$

The test statistic therefore has a value of $K = 14\cdot `r nominator`/`r denominator` = `r 14 * nominator/denominator`$. Is is large enough so that we can claim to have found differences between the groups? You know by heart now that we need to compare this value to a suitable null distribution. The case for the Kruskal Wallis null distribution is similar to what we heard before, i.e. we have two options

1.we enumerate all the possible assignments of the ranks $1-n$ to the groups keeping the group--specific sample sizes $n_i$ the same (for our example this would mean finding all the possibilities to assign the values $1, \ldots , 15$ to three groups each having a size of 5), calculate $K$ each time and see how extreme our observed $K$ is, or,

2. we use an approximation that works best for large sample sizes. Here the appropriate approximate null distribution is the $\chi^2$--distribution with $g-1$ degrees of freedom.

Most software packages use option number two which for our example yields a $p$--value of $`r myKW$p.value`$^[The 'coin' package in R also offers option number 1]. How does that compare to our ANOVA analysis? Recall that we observed a $p$--value of $`r ANOVA.p`$.







## Fisher's exact test

```{r echo = FALSE}

k <- 7
x <- 0:(k)
m <- 10
n <- 10
obs <- 5
FisherProb <- dhyper(x, m, n, k)
options(digits = 4)
```

There is one more non--parametric test on our list, namely Fisher's exact test. From the overview in Table _np_table_ we see that it is an alternative to the $\chi^2$--test. As mentioned above, both the $\chi^2$--test and Fisher's exact test are non--parametric. In fact, they cover the two options that we talked about in the context of the other non--parametric tests in that the $\chi^2$-test uses an approximate null distribution that is inherently based on the Central Limit Theorem (and therefore requires a largish sample size) whereas the procedure that was developed by Sir Ronald Fisher is based on an exact null distribution that can be derived by enumerating all possible outcomes of the test statistic and counting how often they are more extreme than the value we observe in our data.

Let's return to our Martians for the last time. We'll use both sites this time, so that we have 20 individuals. Recall that we put this analysis on hold when we covered the $\chi^2$--test (cf. page _chistoosmaln_) because the smallest group should have no less then 5 observations. But in the context of Fisher's exact test, we can work with this data:

<!-- Observed table of counts}: -->

<!-- Nose&No&Yes&Total\tabularnewline -->
<!-- Red & $2$ & $8$ & $10$\tabularnewline -->
<!-- Blue & $5$ & $5$ & $10$\tabularnewline -->

<!-- Total & $7$ & $13$ & $20$\tabularnewline -->
_FisherExact_
Observed contingency table for Martian colour and presence/absence of antennae in our entire Martian sample - 20 observations in total, taken from both sites.

<!-- % In total, we observed only 20 individuals, and the values in each cell are pretty small, one has only 2 observations, which would be a hint that the $\chi^2$--test isn't the appropriate option. Instead, we can use Fisher's exact test. To do so, we first need to know how likely it is to observe a given contingency table under the null hypothesis of no association between the two variables A and B. Luckily, this is a very old statistical problem that is used a lot in class and usually bores people out of their minds, so we'll be quick! -->

Let's re-label the contingency table, keeping the numbers the same:

<!-- \multicolumn{3}{l}{\emph{Relabelled table}:} \tabularnewline -->

<!-- Colour&Yes&No&Total\tabularnewline -->
<!-- \hline -->
<!-- Red & $2$ & $8$ & $10$\tabularnewline -->
<!-- White & $5$ & $5$ & $10$\tabularnewline -->
<!-- \hline -->
<!-- Total & $7$ & $13$ & $20$\tabularnewline -->


Now, the table could be interpreted as the outcome of an experiment where there are 20 marbles in an urn (yaaawn...), 10 out of these are red and 10 are white. The experimental protocol tells us to draw 7 marbles and the outcome $X$ we're interested in is how many of them are white.^[Note that we can fill in the rest of the table once we know the number of white marbles we draw because the row and column totals are fixed in this experiment!] The observed outcome in our example $X=5$. What's the probability?!

This is a very old statistical problem, i.e. repeated drawing without replacement (because we don't put the marbles we have drawn back into the urn). People have thought about it before so there is a formula we can apply to calculate how likely it is to draw exactly 5 white marbles in 7 draws. The formula is based on the hypergeometric distribution and the calculated probabilities for this example are plotted in Figure _HyperProb_. We see that the resulting probability for $X=5$ is $p_5 = `r FisherProb[6]`$.

```{r echo = FALSE}

fish.df <- data.frame(x = x,
                      prob = dhyper(x, m, n, k))

fish.df$col = MidYellow
# Two sided colours:
# fish.df$col[fish.df$prob <= dhyper(obs, m, n, k)] <- MidRed

# One sided colours:
fish.df$col[fish.df$x >= obs] <- MidRed

fish.df$col <- factor(fish.df$col)

# Draw basic yellow probability plot:
ggplot(fish.df, aes(x = x, y = prob, fill = col)) +
  geom_bar(stat = "identity") +
  geom_vline(xintercept = obs) +
  geom_hline(yintercept = dhyper(obs, m, n, k)) +
  scale_fill_identity(guide = FALSE) +

  scale_x_continuous("Number of white balls",
                     limit = c(-0.5, 7.5),
                     labels = seq(0,7,1),
                     breaks = seq(0,7,1),
                     expand = c(0,0)) +
  scale_y_continuous("Probability",
                     limit = c(0, 0.35),
                     labels = seq(0, 0.35, 0.05),
                     breaks = seq(0, 0.35, 0.05),
                     expand = c(0,0))  +
  theme(panel.background = element_blank(),
        plot.background = element_blank(),
        legend.background = element_blank(),
        legend.key = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.ticks.y = element_line(colour = "black"),
        axis.ticks.x = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
  )
```
_HyperProb_
The hypergeometric distribution. Drawing 7 balls from an urn with 10 red and 10 white balls. The $x$--axis refers to the number of white balls, the $y$--axis shows the corresponding point probability. The red bars correspond to the probability of drawing 5 or more white balls which is the (one--sided) $p$--value.

Note that in the context of the urn analogy, the inherent assumption is that we are not influenced by the colour of the marbles when we draw them. This corresponds to the null hypothesis of Fisher's exact test. Just like in the $\chi^2$--test, it is that both variables ('drawn (yes/no)' and 'colour (yes/no)') are not associated. So Figure _FisherExact_ actually shows us the null distribution of our experiment! Also we already know that the obseved value of the test statistic is $X=5$.

Now we are only one step away from the $p$--value - the probability to observe data that differ from the null hypothesis _at least_ as much as the ones at hand, assuming $H_0$ is true. We need to come up with the other possible contingency tables that differ even more from the null hypothesis of no association than the one above. In our example, this is quite easy and the only other two options are shown in Figure _FisherExact_. The reason there are only two more options is that we don't want to change the experimental setup. i.e. in the urn there are 10 red and 10 white marbles and we draw 7 times. The originally observed number of 5 white marbles is quite extreme already, only 6/7 and 7/7 would be more extreme. The corresponding probabilities are $p_6 = `r FisherProb[7]`$ and  $p_7 = `r FisherProb[8]`$. You can also see these in Figure _FisherExact_.


<!-- More extreme cases: -->

<!-- Colour&Yes&No&Total\tabularnewline -->
<!-- \hline -->
<!-- Red & $1$ & $9$ & $10$\tabularnewline -->
<!-- White & $6$ & $4$ & $10$\tabularnewline -->
<!-- \hline -->
<!-- Total & $7$ & $13$ & $20$\tabularnewline -->


<!-- Colour&Yes&No&Total\tabularnewline -->
<!-- \hline -->
<!-- Red & $0$ & $10$ & $10$\tabularnewline -->
<!-- White & $7$ & $3$ & $10$\tabularnewline -->
<!-- \hline -->
<!-- Total & $7$ & $13$ & $20$\tabularnewline -->

_FisherExact_
Contingency tables that are more extreme (i.e. more deviant from the null hypothesis of no association) than the table we observed.




This is all we need to calculate the $p$--value of Fisher's exact test which for our example is $p_5 + p_6 + p_7 = `r sum(dhyper(obs:k, m, n, k))`$. Therefore, for our example we cannot reject the $H_0$ of no association at the 5\% significance level. Note that this result also applies for the Martian data (all we did to transform it into the urn context was relabel the variable, that doesn't change the result). Also note that we considered a one--sided test, if we wanted to account for a two--sided test we would have doubled our $p$--value.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

<!-- % \section{Exact versus Approximation} -->
<!-- %  -->
<!-- %     There are two things worth noting here. First, the Wilcoxon signed--rank test assumes that the observed differences are symmetrically distributed around the (population) median. So, we are not entirely free of assumptions.  -->
<!-- %      -->
<!-- %      -->
<!-- % There are two methods for deriving the null distribution for our test statistic. The first method is _exact_. We used this above when we calculated _all_ possible values that the test statistic could take and from that derived a null distributuion and $p$--values. There was no approximation. The second method uses an _approximate null distribution_, which is based on the central limit theorem. These two methods will make a reappearance for all the tests that we will see here. -->
<!-- % Exact tests seem like a good thing, but for large samples, finding _all_ possibilities can take a while - even for computers, and if $n$ is large, the approximation method yields very similar results a lot quicker. Some statistical software packages offer both methods, others will just use one -- maybe without telling you. If you want to know what's going on, it's worth reading the help files of your particular package! _approx_NB_ -->


## Some final words

To summarise, non--parametric tests have many advantages over parametric tests. They can be used with small sample sizes and don't require Normally-distributed data, i.e. they are able to work with "dirty" data, which may contain outliers and other difficult problems. However, because they also make fewer assumptions about the nature of the data, they operate on a coarser scale, making them less sensitive for detecting interesting patterns in the data. Everything in life is a trade--off, and Statistics is no exception!

Most non--parametric tests have either an exact or an approximate version that yield very similar results if the sample size is large. If not, the exact version is preferable. Unfortunately, this is the situation where the increased coarseness of the non--parametric tests really makes a differences. As we saw, sometimes the $p$--values that can be obtained in a non--parametric test comes in very big steps, and sometimes it can be impossible to reject the null hypothesis at the 5\% significance level regardless of the data, just because the lowest obtainable $p$--value with a given sample size is above 0.05!

While this may seem annoying, it's quite honest. For small sample sizes, the evidence contained in your dataset is sparse, and unless you are willing to make strong assumptions that justify parametric methods (where small $p$--values can be potentially obtained even with small $n$), there is just not much you can say about the world if $n$ is small. Rejecting a null hypothesis with only a 5\% probability of being wrong is sometimes just not possible if you don't have enough information.
The following quote sums up this point here quite nicely

> Finally, researchers should always judge the credibility of their findings and should remember that extraordinary claims require extraordinary evidence. ^[Practical Assessment, Research \& Evaluation, Vol 18, No 10 Page 2 De Winter; t-test with extremely small Ns.]