Chapter 7 Hypothesis testing

7.1 Refresher

Copy link

We have flipped a coing 100 times and have obtained the following results.

## flips
##  H  T 
## 62 38

If the coin was unbiased we expect roughly 50 heads. Can we make any claims regarding the biaised or unbiaised nature of that coin?

A coin toss can be modeled by a binomial distribution. The histogram below shows the binomial statistic for 0 to 100 heads; it represents the binomial density of an unbiased coin. The vertical line shows the number of head observed.

Figure 7.1: Binomial density of an unbiased coin to get 0 to 100 heads. The full area of the histogram sums to 1.

Above, we see that the most likely outcome would be 50 heads with a probability of 0.0795892. No head and 100 heads have a probability of \(7.8886091\times 10^{-31}\).

We set

• \(H_0\): the coin is fair
• \(H_1\): the coin is baised

If 62 isn’t deemed too extreme, then we won’t reject \(H_0\) and conclude that the coin in fair. If 62 is deemed too extreme, then we reject \(H_0\) and accept \(H_1\), and conclude that the coin is baised.

To define extreme, we set \(\alpha = 0.05\) and rejet \(H_0\) if our result is outside of the 95% most probably values.

Figure 7.2: Binomial density of for an unbiased coin to get 0 to 100 heads. The areas in red sum to 0.05.

We can also compute the p-value, the tells us how likely we are to see such a extreme or more extreme value under \(H_0\).

binom.test(x = 62, n = 100, p = 0.5)

## 
##  Exact binomial test
## 
## data:  62 and 100
## number of successes = 62, number of trials = 100, p-value = 0.02098
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.5174607 0.7152325
## sample estimates:
## probability of success 
##                   0.62

Whenever we make such a decision, we will be in one of the following situations:

See below for a step by step guide to this example.

7.2 A biological example

Copy link

A typical biological example consists in measuring a gene of interest in two populations of interest (say control and group), represented by biological replicates. The figure below represents the distribution of the gene of interest in the two populations, with the expression intensities of triplicates in each population.

Figure 7.3: Expression of a gene in two populations with randomly chosen triplicates.

We don’t have access to the whole population and thus use a sample thereof (the replicated measurements) to estimate the population parameters. We have

We set our hypothesis as

• \(H_0\): the means of the two groups are the same, \(\mu_{1} = \mu_{2}\).
• \(H_1\): the means of the two groups are different, \(\mu_{1} \neq \mu_{2}\).

and calculate a two-sided, two-sample t-test (assuming unequal variances) with

\[ t = \frac{ \bar X_{1} - \bar X_{2} } { \sqrt{ \frac{ s_{1}^{2} }{ N_{1} } + \frac{ s_{2}^{2} }{ N_{2} } } } \]

where \(\bar X_i\), \(s_{i}^{2}\) and \(N_i\) are the mean, variance and size of samples 1 and 2.

A t-test has the following assumptions:

• the data are normally distributed;
• the data are independent and identically distributed;
• equal or un-equal (Welch test) variances.

Note that the t-test is robust to deviations88 All models are wrong. Some are useful..

In R, we do this with

s1

## [1] 5.327429 5.723787 7.870450

s2

## [1] 11.10576 11.19393 13.57260

t.test(s1, s2)

## 
##  Welch Two Sample t-test
## 
## data:  s1 and s2
## t = -5.0003, df = 3.998, p-value = 0.007499
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -8.788162 -2.512255
## sample estimates:
## mean of x mean of y 
##  6.307222 11.957431

In practice, we would apply moderated versions of the tests, such as the one provided in the limma package and also widely applied to RNA-Seq count data.

7.3 A more realistic biological example

Copy link

Let’s now use the tdata1 dataset from the rWSBIM1322 package that provide gene expression data for 100 genes and 6 samples, three in group A and 3 in group B.

library("rWSBIM1322")
data(tdata1)
head(tdata1)

##                 A1         A2         A3         B1         B2        B3
## feature1  429.6721   433.3133  1806.2773  3534.6396  4541.2757 1139.8326
## feature2 1320.8827  1085.5660 10643.9541   239.7528 14177.4412  678.5271
## feature3  321.5151   289.6480  9236.7965 15745.7491   443.8737  198.5716
## feature4 9348.7749  1274.8021   647.2523   601.6313   595.6266 1040.1568
## feature5 1616.9009   413.2385    43.0974 10143.6595   575.0549 4048.5851
## feature6  328.3438 11866.2342 32661.9472  8331.8745   608.3194 9331.7692

log_tdata1 <- log2(tdata1)
par(mfrow = c(1, 2))
limma::plotDensities(tdata1)
limma::plotDensities(log_tdata1)

We are now going to apply a t-test to feature (row) 73, comparing the expression intensities in groups A and B. As we have seen, this can be done with the t.test function:

x <- log_tdata1[73, ]
t.test(x[1:3], x[4:6])

## 
##  Welch Two Sample t-test
## 
## data:  x[1:3] and x[4:6]
## t = 0.1728, df = 3.999, p-value = 0.8712
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.444807  3.902010
## sample estimates:
## mean of x mean of y 
##  10.14391   9.91531

my_t_test <- function(x) {
    t.test(x[1:3], x[4:6])$p.value
}
pvals <- apply(log_tdata1, 1, my_t_test)
table(pvals < 0.05)

## 
## FALSE  TRUE 
##    98     2

head(sort(pvals))

##  feature19  feature61  feature20  feature84   feature1   feature7 
## 0.02967230 0.03595002 0.08576252 0.08958506 0.10645476 0.10832969

To answer these questions, let’s refer to this xkcd cartoon that depicts scientists testing whether eating jelly beans causes acne.

Figure 7.4: Do jelly beans cause acne? Scientists investigate. From xkcd.

The data that was used to calculate the p-values was all drawn from the same distribution \(N(10, 2)\).

As a result, we should not expect to find a statistically significant result, unless we repeat the test enought times. Enough here depends on the \(\alpha\) we set to control the type I error. If we set \(\alpha\) to 0.05, we accept that rejecting \(H_{0}\) in 5% of the extreme cases where we shouldn’t reject it. This is an acceptable threshold that however doesn’t hold when we repeat the test many time.

An important visualistion when performing statistical test repeatedly, is to visualise the distribution of computed p-values. Below, we see the histogram of the tdata1 data and 100 values drawn from a uniform distribution between 0 and 1. Both are very similar; they are flat.

par(mfrow = c(1, 2))
hist(pvals)
hist(runif(100))

Figure 7.5: Distribution of p-values for the tdata1 dataset (left) and 100 (p-)values uniformely distributed between 0 and 1 (right).

Below we see the expected trends of p-values for different scenarios. This example comes from the Variance explained.

Figure 7.6: Expected trends of p-values for different scenarios (source).

In an experiment with enough truly differentially expression features, on expects to observe a substantial increase of small p-values (anti-conservative). In other words, we expect to see more small p-values that at random, or when no statistically significant are present (uniform). All other scenarios warrant further inspection, as they might point to issue with the data of the tests.

When many tests are preformed, the p-values need to be adjusted, to take into account that many tests have been performed.

7.4 Adjustment for multiple testing

Copy link

There are two classes of multiple testing adjustment methods:

• Family-wise error rate (FWER) that gives the probability of one or more false positives. The Bonferroni correction for m tests multiplies each p-value by m. One then checks if any results still remains below significance threshold.

• False discovery rate (FDR) that computes the expected fraction of false positives among all discoveries. It allows us to choose n results with a given FDR. Widely used examples are Benjamini-Hochberg or q-values.

The figure below illustrates the principle behind the FDR adjustment. The procedure estimate the proportion of hypothesis that are null and then adjust the p-values accordingly.

Figure 7.7: Principle behind the false discovery rate p-value adjustment (source).

Back to our tdata1 example, we need to take this into account and adjust the p-values for multiple testing. Below we apply the Benjamini-Hochberg FDR procedure using the p.adjust function and confirm that none of the features are differentially expressed.

adj_pvals <- p.adjust(pvals, method = "BH")

any(adj_pvals < 0.05)

## [1] FALSE

min(adj_pvals)

## [1] 0.899736

7.5 Result visualisation

Copy link

library("rWSBIM1322")
data(tdata4)

my_t_test <- function(x)
    t.test(x[1:3], x[4:6])$p.value

pvals <- apply(tdata4, 1, my_t_test)
hist(pvals)

table(pvals < 0.05)

## 
## FALSE  TRUE 
##    91     9

adjp <- p.adjust(pvals, method = "BH")
head(sort(adjp))

##  feature87  feature34  feature39  feature68  feature43  feature51 
## 0.02927868 0.09032688 0.09032688 0.09032688 0.21467090 0.25799641

my_lfc <- function(x)
    mean(x[1:3]) - mean(x[4:6])
lfc <- apply(tdata4, 1, my_lfc)
hist(lfc)

We generally want to consider both the statistical significance (the p-value) and the magnitude of the difference of expression (the fold-change) to provide a biological interpretation of the data. For this, we use a volcano plot as shown below. The most interesting features are those towards to top corners given that they have a small p-value (i.e. a large value for -log10(p-value)) and large (in absolute value) fold-changes.

Figure 7.8: A volcano plot.

plot(lfc, -log10(adjp))
grid()
abline(h = -log10(0.05))
abline(v = c(-1, 1))

Chapter 6 on testing of Modern Statistics for Modern Biology (⊕Holmes and Huber 2019Holmes, Susan, and Wolfgang Huber. 2019. Modern Statistics for Modern Biology. Cambridge Univeristy Press.), and in particlar section 6.5, provides additional details about the t-test.

The t-test comes in multiple flavors, all of which can be chosen through parameters of the t.test function. What we did above is called a two-sided two-sample unpaired test with unequal variance. Two-sided refers to the fact that we were open to reject the null hypothesis if the expression in the first group was either larger or smaller than that in the second one.

Two-sample indicates that we compared the means of two groups to each other; another option is to compare the mean of one group against a given, fixed number.

Unpaired means that there was no direct 1:1 mapping between the measurements in the two groups. If, on the other hand, the data had been measured on the same patients before and after treatment, then a paired test would be more appropriate, as it looks at the change of expression within each patient, rather than their absolute expression.

Equal variance refers to the way the statistic is calculated. That expression is most appropriate if the variances within each group are about the same. If they are very different, an alternative form, called the Welch t-test exist.

An important assumption when performing a t-test is the independence assumption among the observations in a sample. An additional exercise below demonstrate the risk of dependence between samples.

Finally, it is important to highlight that all the methods and concepts described above will only be relevant if the problem is properly stated. Hence the importance of clearly laying out the experimental design and stating the biological question of interest before running experiment. Indeed, quoting the mathematician Richard Hamming:

It is better to solve the right problem the wrong way than to solve the wrong problem the right way.

To this effect, a type III error has been defined as errors that occur when researchers provide the right answer to the wrong question.

7.6 A word of caution

Copy link

Summary of statistical inference:

1. Set up a model of reality: null hypothesis \(H_0\) (no difference) and alternative hypothesis \(H_1\) (there is a difference).

2. Do an experiment to collect data.

3. Perform the statistical inference with the collected data.

4. Make a decision: reject \(H_0\) if the computed probability is deemed too small.

What not to do:

7.7 Additional exercises

Copy link

set.seed(2)
n <- 100
p <- 0.59
flips <- sample(c("H", "T"), size = n,
                replace = TRUE,
                prob = c(p, 1 - p))

table(flips)

## flips
##  H  T 
## 62 38

library("dplyr")
num_heads <- sum(flips == "H")
binomial_dens <-
    tibble(k = 0:n) %>%
    mutate(p = dbinom(k, size = n, prob = 0.5))

library("ggplot2")
ggplot(binomial_dens, aes(x = k, y = p)) +
    geom_bar(stat = "identity") +
    geom_vline(xintercept = num_heads)

alpha <- 0.05

binomial_dens %>%
    filter(p == max(p))

## # A tibble: 1 × 2
##       k      p
##   <int>  <dbl>
## 1    50 0.0796

binomial_dens %>%
    filter(p == min(p))

## # A tibble: 2 × 2
##       k        p
##   <int>    <dbl>
## 1     0 7.89e-31
## 2   100 7.89e-31

binomial_dens <-
    arrange(binomial_dens, p) %>%
    mutate(reject = (cumsum(p) <= alpha))

ggplot(binomial_dens) +
    geom_bar(aes(x = k, y = p, fill = reject), stat = "identity") +
    scale_fill_manual(
        values = c(`TRUE` = "red", `FALSE` = "#5a5a5a")) +
    geom_vline(xintercept = num_heads, col = "blue") +
    theme(legend.position = "none")

binom.test(x = num_heads, n = n, p = 0.5)

## 
##  Exact binomial test
## 
## data:  num_heads and n
## number of successes = 62, number of trials = 100, p-value = 0.02098
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.5174607 0.7152325
## sample estimates:
## probability of success 
##                   0.62