---
title: Simple hypothesis testing in R
author: Aaron Lun, Catalina Vallejos
date: 11 November 2016
---

# Testing ratios in contigency tables

## Using Fisher's exact test

We'll be using the classic tea-tasting experiment. 
From `?fisher.test`:

> A British woman claimed to be able to distinguish whether milk or
> tea was added to the cup first.  To test, she was given 8 cups of
> tea, in four of which milk was added first.  The null hypothesis
> is that there is no association between the true order of pouring
> and the woman's guess, the alternative that there is a positive
> association (that the odds ratio is greater than 1).

Let's set up the data set first.
Below is data for 8 cups:

```{r}
true.first <- rep(c("Milk", "Tea"), each=4)
guess.first <- rep(c("Milk", "Tea", "Milk", "Tea"), c(3,1,1,3))
```

Can you summarize this into a 2-by-2 contingency table?

```{r}
tea.data <- table(true.first, guess.first) 
tea.data 
```

For 2-by-2 tables, Fisher's exact test works by testing whether the "odds ratio" is equal to 1.
The "odds ratio" is calculated as a ratio of odds (obviously) between our two experimental conditions, i.e., milk or tea added first.
What are the odds that the guess was milk in the cases where milk was added first? 3:1.
What are the odds that the guess was milk in the cases where tea was added first? 1:3.
Thus, the sample odds ratio between milk/tea first is 3:1/1:3 = 9.

Now, apply `fisher.test` to the table.

```{r}
f.out <- fisher.test(tea.data) 
f.out 
```

How do I extract the _p_-value?
How do I extract the odds ratio? (__Note:__ this is slightly different from the sample odds ratio.)
Does it matter if I switch the rows and columns?

```{r}
f.out$p.value # To get the p-value 
f.out$estimate # To get the odds ratio 

# We can use t() to transpose the table: 
fisher.test(t(tea.data))  
```

<!-- 
Transposition has no effect on the test, because the nature of the null hypothesis is unchanged. 
The odds ratio is the same regardless of whether it is computed using odds per column or per row.
--> 

A key assumption is:

- Observations are **independent**, randomly sampled from a population

Here, the "population" refers to the set of all possible observations taken from this particular woman (not from the population of all women!). 
Each observation (i.e., the guess for each cup) is assumed to be independent of the others.

_Note:_ What is the null hypothesis for larger contingency tables?
That the experimental condition does not affect the guesses, i.e., that the distribution of counts across columns is the same for all rows. 

```{r}
true.first <- rep(c("Milk", "Tea", "Sugar"), each=6)
guess.first <- sample(c("Milk", "Tea", "Sugar"), 18, replace=TRUE)
tea.data.exp <- table(true.first, guess.first)
fisher.test(tea.data.exp) # What is this testing?
```

## Using Pearson's Chi-squared test of independence

Fisher's exact test computes the p-value exactly for small sample sizes.
For large sample sizes, we can use Pearson's Chi-squared test instead.
Let's mock up some data by imagining we repeated our tea experiment 100 times:

```{r}
true.first <- rep(c("Milk", "Tea"), each=4*100)
guess.first <- sample(c("Milk", "Tea", "Milk", "Tea"), 800, replace=TRUE)
tea.data.more <- table(true.first, guess.first)
tea.data.more
```

Try applying Fisher's exact test and Pearson's Chi-squared test on this data.
Do the p-values match up between the tests?

```{r}
fisher.test(tea.data.more) 
chisq.test(tea.data.more) 
```

In addition to the **independence** assumption (as in the Fisher's exact test), we require that:

- at least 80% of the cells have an **expected** frequency of 5 or greater

- none of the cells have an **expected** frequency less than 1 

This is necessary for the distribution of the test statistic to be accurately approximated with a chi-squared distribution.
What happens if there aren't enough observations (try it out on `tea.data`)?

```{r}
fisher.test(tea.data) 
chisq.test(tea.data) 
```

What happens if we have really large tables?

```{r}
true.first <- rep(c("Milk", "Tea", "Sugar"), each=600)
guess.first <- sample(c("Milk", "Tea", "Sugar"), 1800, replace=TRUE)
tea.data.exp.more <- table(true.first, guess.first)
try(fisher.test(tea.data.exp.more)) 
chisq.test(tea.data.exp.more) 
```

The Chi-squared test is also useful for testing goodness of fit.
Say we have a bunch of F2 mice.
Do their genotypes follow the 1:2:1 Mendelian ratio (hint: use the `p` argument in `chisq.test`)?

```{r}
genotypes <- rep(c("AA", "Aa", "aa"), c(500, 600, 100))
# Create a table of genotypes:
genotab <- table(genotypes) 
# Create a vector of relative expected frequencies:
genoexp <- c(0.25, 0.5, 0.25)
chisq.test(genotab, p=genoexp) 
```

# Using _t_-tests of various flavours

## One sample t-test

This is useful for comparing our observations against an expected value under the null hypothesis.
For example, I set the thermostat for my office to 25 degrees Celsius.
To check whether it's at the intended temperature, I measure the temperature every day for a month.

```{r}
temp <- c(22.8, 22.7, 25.9, 22.0, 21.2, 24.4, 22.1, 22.6, 23.3, 22.2, 
          22.5, 24.1, 22.4, 25.1, 24.1, 23.3, 23.3, 22.5, 23.8, 22.4, 
          22.5, 23.4, 23.5, 22.5, 25.8, 23.8, 22.0, 22.0, 23.2, 22.0)
```

My null hypothesis is that the temperature is indeed 25 degrees Celsius.
My alternative hypothesis is that, well, it's not.
How can I apply a one-sample t-test to this data (hint: look at `mu` in `t.test`)?

```{r}
t.out <- t.test(temp, mu=25) 
t.out 
```

How can I extract the _p_-value, or the _t_-statistic?

```{r}
t.out$p.value 
t.out$statistic 
```

The one-sample $t$-test assumes that:

- Observations are **independent**, randomly sampled from a population

- Observations are **normally distributed**

When the sample size is large (rule of thumb $n > 30$), the Central Limit Theorem allows us to use this test (even when the data is not normal!)

Let's say I don't want to waste energy during winter, so I want to check whether the temperature is above 25 degrees Celsius.
Here, my null hypothesis is that the temperature is below or equal to 25 degrees.
My alternative hypothesis is that the temperature is above 25 degrees.
This is a one-sided, one-sample t-test (whereas the previous test above was a two-sided, one-sample t-test).
How can I do this with `t.test`?

```{r}
t.test(temp, mu=25, alternative="greater") 
```

## Two sample t-test

This is useful for comparing between observations from two groups.
For example, let's say I run a chicken farm.
I randomly select 20 chickens and feed them with soybean.
I randomly select another 30 chickens and feed them with sunflower seeds.
I measure their weights after 6 weeks:

```{r}
on.soy <- c(188, 229, 192, 207, 172, 151, 188, 173, 209, 158, 
            211, 201, 205, 249, 214, 171, 198, 192, 161, 179)
on.sun <- c(213, 262, 237, 223, 223, 268, 216, 276, 258, 269,
            227, 252, 243, 286, 247, 277, 242, 233, 228, 211, 
            230, 254, 268, 242, 223, 242, 224, 279, 227, 251)
```

I want to see if there is a difference in weight between the feed types.
My null hypothesis is that the average weights are equal between the feed types.
My alternative hypothesis is that they are not.
How can I apply a two-sample t-test to this data?

```{r}
t.test(on.sun, on.soy) 
```

The two-sample $t$-test assumes that

- **Within** each group, observations are **independent** and randomly drawn from a population

- Observations are independent **between** the groups

- Observations are normally distributed **within** each group

As in the one-sample t-test, the Central Limit Theorem allows us to use this test for non-normal data when the sample size is large (rule of thumb: > 30 observations per group).

The groups don't have to be of equal size (i.e., "unbalanced"), though greater power is usually achieved if they are equal.

__Note:__ By default, `t.test` uses Welch's t-test, which allows the variances in each of the two groups to be unequal.
If `var.equal=TRUE`, the code will use Student's t-test, which makes the additional assumption of equal variances in each of the two groups.
It's usually best to stick to the default, which is more robust.

```{r}
t.test(on.sun, on.soy, var.equal=TRUE) 
```

Let's say I currently feed all my chickens on soybean.
I want to know whether sunflower feed provides more weight gain - I don't care to know if it results in less weight. 
Here, my null hypothesis is that the weight with sunflower feed is less than or equal to the weight with soybean.
My alternative hypothesis is that the weight with sunflower feed is greater than that with soybean.
This is a one-sided, two-sample t-test (whereas the previous test was a two-sided, two-sample t-test).
How can I do this with `t.test`?

```{r}
t.test(on.sun, on.soy, alternative="greater") 
# Note that 'greater' checks for "the first argument ('on.sun') > second ('on.soy')". 
```

## Paired samples t-test

This is useful for comparing between groups where the samples are paired in some manner.
For example, let's say I randomly sample 15 patients and collect healthy and diseased tissue from each.
I measure the expression of a gene in each of healthy/disease samples (e.g., via qPCR):

```{r}
                    #1   #2   #3   #4   #5   #6   #7    #8   #9   #10  #11  #12  #13  #14  #15
healthy.exprs <- c(6.9, 8.2, 9.0, 7.7, 6.5, 7.7, 5.5,  7.7, 7.6,  8.2, 6.5, 6.3, 8.7, 6.2, 8.4)
disease.exprs <- c(8.3, 9.5, 9.8, 9.9, 6.6, 8.1, 9.5, 10.0, 8.7, 10.0, 7.5, 8.0, 8.5, 5.0, 8.5)
```

I want to know if the gene is differentially expressed upon disease.
My null hypothesis is that the average expression of the gene is the same between disease and healthy groups.
Could I use the two-sample t-test?
The observations are not independent between groups, because each pair of healthy/diseased tissue comes from the same patients; so no. 
How could I do a paired-sample t-test?

```{r}
t.test(healthy.exprs, disease.exprs, paired=TRUE) 
```

Some of the assumptions are simular to the ones for two-sample $t$-test, namely:

- **Within** each group, observations are **independent** and randomly drawn from a population

However, the paired $t$-test does not assume independence between the groups. Instead, it assumes that:

- Pairs of samples are randomly drawn from a population (i.e. each pair is **independent** from all the others)

- Differences between paired observations are normally distributed.

As before, the Central Limit Theorem allows us to use this test when the sample size is large (rule of thumb: at least 30 pairs of observations).

Just like the other t-tests, this can also be two-sided or one-sided, depending on what alternative hypothesis is of interest.

## Are the data normal?

A key assumption of all the previous tests was that that data were normally distributed.
We can check this with histograms:

```{r}
data1 <- rnorm(1000)
data2 <- runif(1000)
par(mfrow=c(1,2)) 
hist(data1) # looks normal
hist(data2) # doesn't look normal.
```

Alternatively, rather than trying to match the shape of the histogram to the shape of the normal distribution, we can use `qqnorm` and `qqline`.
These plot the quantiles of the observed distribution (i.e., the value at 10%, 25%, 50%, etc.) against the expected quantiles for a normal distribution.
A normally-distributed set of observations should lie on the diagonal of the quantile-quantile plot. 

```{r}
par(mfrow=c(1,2)) 
qqnorm(data1) 
qqline(data1) 
qqnorm(data2) 
qqline(data2) 
```

How do I do this if I have two groups?

```{r}
g1 <- rnorm(1000, mean=1)
g2 <- rnorm(1000, mean=2)
par(mfrow=c(1,2)) 
qqnorm(g1) 
qqline(g1) 
qqnorm(g2) 
qqline(g2) 
```

What if the data was paired?

```{r}
disease <- rnorm(1000, mean=1)
healthy <- rnorm(1000, mean=2)
effect <- disease - healthy 
qqnorm(effect)
qqline(effect) 
```

Major deviations from the diagonal indicate there may be problems with using t-tests.
(Though this is generally difficult to check with few observations in most experimental datasets.
Formal tests for deviation from normality exist, but they are less useful than one would expect.)
Whenever possible, it's generally safer to use non-parametric tests.

# Non-parametric tests of various flavours

## Using the Mann-Whitney U (two-sample Wilcoxon) test

Let's re-use our chicken feeding example:

```{r}
on.soy <- c(188, 229, 192, 207, 172, 151, 188, 173, 209, 158, 
            211, 201, 205, 249, 214, 171, 198, 192, 161, 179)
on.sun <- c(213, 262, 237, 223, 223, 268, 216, 276, 258, 269,
            227, 252, 243, 286, 247, 277, 242, 233, 228, 211, 
            230, 254, 268, 242, 223, 242, 224, 279, 227, 251)
```

Our null hypothesis is that these observations come from the same distribution.
By implication, this means that there is no shift in the average weight between feed types.
We can test this hypothesis using `wilcox.test`:

```{r}
wilcox.test(on.soy, on.sun) 
```

Key assumptions include:

- observations are randomly sampled from the population corresponding to each group
- observations are independent within and between groups

No normality required.

__Note:__ The test outcome can be interpreted as a shift in medians, but only if the shape of the distribution does not change between groups.
It's worthwhile looking at the data as a sanity check during interpretation:

```{r}
boxplot(list(Soybean=on.soy, Sunflower=on.sun), ylab="Weight (g)") 
```

(It also is possible to do one-sided tests here, but they are somewhat difficult to interpret if the distribution changes.)

## Using the Wilcoxon signed-rank test

Let's re-use our patient example:

```{r}
                    #1   #2   #3   #4   #5   #6   #7    #8   #9   #10  #11  #12  #13  #14  #15
healthy.exprs <- c(6.9, 8.2, 9.0, 7.7, 6.5, 7.7, 5.5,  7.7, 7.6,  8.2, 6.5, 6.3, 8.7, 6.2, 8.4)
disease.exprs <- c(8.3, 9.5, 9.8, 9.9, 6.6, 8.1, 9.5, 10.0, 8.7, 10.0, 7.5, 8.0, 8.5, 5.0, 8.5)
```

Our null hypothesis is that the difference in expression is symmetrically distributed around zero.
This indicates that there is no systematic difference between healthy and diseased tissue. 
To test this, we use `wilcox.test` but with `paired=TRUE`:

```{r}
wilcox.test(healthy.exprs, disease.exprs, paired=TRUE) 
```

Key assumptions include:

- pairs are randomly sampled from the population
- each pair is independent of the others

Again, no normality required.

__Note:__ Asymmetry can cause rejection without any actual shift in the median or mean of the differences.
Thus, it's useful to examine the distribution of differences to check for a shift:

```{r}
boxplot(disease.exprs - healthy.exprs, ylab="Difference in expression") 
```

# Session information

We save the session information into the report for posterity.

```{r}
sessionInfo()
```
