This document is released with the following license: CC-BY 4.0

How to use this document

This document is meant to be read and played with interactively. Whenever you see the line ## Parameters in the code, you are encouraged to try different values for the parameters, to see how the results change.

What is inferential statistics?

Inferential statistics is the process of using sample data to draw conclusions and make generalizations about a larger population while accounting for uncertainty and random variation.

In contrast to descriptive statistics, we are trying to go beyond the information presented by the sample, and try to answer hypotheses about our population.

One of the most important tools of inferential statistics is hypothesis testing. In its simplest and most common form, we want to use a statistical method to make a decision about a hypothesis.

However, more common statistical methods in our field (frequentist approach) do not assign probability to the validity of a hypothesis, but rather focus on calculating the probability (p-value) that our data points can be explained by a certain assumed distribution (null hypothesis).

What is a p-value?

The p-value is the probability of observing your data (or more extreme data) assuming the null hypothesis is true.

A p-value can be calculated with any distribution, depending on the test statistic.

The following represents how the p-value is calculated based on an assumed distribution. Taken a value (on the x axis), the area of the distribution curve that is further from that specific point is its p-value.

## Parameters

value <- 1.5

##


# Create normal distribution data
x <- seq(-4, 4, length.out = 1000)
y <- dnorm(x)

# Create data frame for plotting
df <- data.frame(x = x, y = y)

p_value <- (1-pnorm(value))

# Create the plot
ggplot(df, aes(x = x, y = y)) +
  geom_line() +
  geom_area(data = subset(df, x >= value), fill = "red", alpha = 0.3) +
  labs(x = "Test Statistic", y = "Density",
       title = paste("Normal Distribution with the p-value for",
                     value, "represented as area (p=", round(p_value, 3), ")")) +
  theme_minimal()

In the example above, we used a normal distribution as the test statistic. As we will see, this is not the only (or even the most common) choice.

The t-test

The t-test is one of the simplest significance tests that we can use. It was developed by William Sealy Gosset while working on quality control at the Guinness Brewery and it addresses the problem of small sample sizes.

A t-test tells us the probability that the mean of a sample comes from a specific distribution (null hypothesis).

Remember! In frequentist statistics, the test gives us the probability of the data under the null hypothesis and not the probability of the null (or alternative) hypothesis!

Our test statistics therefore concerns the distribution of means, not the distribution of the population. From this mean, we can calculate a t-value using the formula:

\[t=\frac{\bar{x}-\mu}{\hat{\sigma} / \sqrt{n}}\] where \(\bar{x}\) is the mean of our sample, \(\mu\) is the mean of the null distribution (which can be 0 for simplicity), \(n\) is the sample size, and \(\hat{\sigma}\) is the estimated standard deviation of the sample: \[\hat{\sigma} = \sqrt{\frac{1}{n-1} \sum_i (x_i - \bar{x})^2}\] The t-value is basically the distance between the mean of our sample from our expected value, expressed in “units” of expected standard error \(\hat\sigma/\sqrt n\)

Why does this work? Why the scaling factor \(\sqrt{n}\)?

This works because we know that, under appropriate conditions, the distribution of the means of a sampling of any other population distribution converges to a normal distribution of mean \(\mu\) and standard deviation \(\sigma/\sqrt{n}\), which is called standard error (central limit theorem).

Important! The distribution of the population doesn’t necessarily need to be normal! Only the distribution of the means need to be.

The central limit theorem

We will now demonstrate that for different population distributions, the means converge to a normal distribution.

## Parameters

# available distributions to test:
# normal_distribution, uniform_distribution, biexp_distribution, skewed_distribution

generate_data <- uniform_distribution
sample_size <- 10  # Size of each sample

##

The t-test procedure

For the simple case of a one-sample t-test:

  1. Measure a sample
  2. Calculate the mean \(\bar{x} = \frac{\sum{x_i}}{n}\) and the estimated standard deviation \(\hat{\sigma} = \sqrt{\frac{1}{n-1} \sum_i (x_i - \bar{x})^2}\)
  3. Calculate the t-value \(t=\frac{\bar{x}-\mu}{\hat{\sigma} / \sqrt{n}}\)
  4. Calculate the p-value by calculating the likelihood of observing the calculated t-value under the assumption of the null hypothesis.

The t-values follow a distribution called the t-distribution. This distribution depends on the sample size, and for small sample sizes it has very “fat” tails. The t-distribution converges to a normal distribution for large sample sizes.

# Create histogram with t-distribution overlay
t_dist_plot <- ggplot(data.frame(t_stats = sample_t_stats), aes(x = t_stats)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "blue", alpha = 0.3) +
  stat_function(fun = dt, 
                args = list(df = sample_size - 1),
                color = "red", size = 1) +
  stat_function(fun = dnorm, 
                color = "darkgreen", size = 1) +
  labs(title = paste("Distribution of T-Statistics (sample size", sample_size, ")"),
       subtitle = "Red line shows theoretical t-distribution, green shows the normal distribution (sample size=∞)") +
  xlim(-6,6) +
  theme_minimal()

t_dist_plot

t-test: summary

A t-test is a test for the mean. We care about the distribution of the means, not necessarily of the population. The means are normally distributed when:

  • the population is normally distributed; or
  • the sample size is large enough (regardless of the population distribution).

The sample size is important. Our t value is scaled by the square root of the sample size. Which means that the same difference in means, with the same standard deviation, will give a larger t-value (and therefore a larger normalized distance from the null hypothesis) with a larger sample size. This concept is related to power and sample size calculations.

Confidence intervals

The t-distribution can be used to calculate confidence intervals for our means. A X% confidence interval comes from a procedure that, if repeated many times with different random samples, would contain the true population mean X% of the time. When talking about means, we can use the t-distribution to build a confidence interval in the following way:

\[CI(X\%) = \bar x \pm t(X\%) \cdot (\hat\sigma/\sqrt n)\] where \(t(X\%)\) is the critical value that corresponds to the desired confidence percentage. In formulas:

## Parameters

population_mean <- 4
population_sd <- 1
sample_size <- 10
confidence <- 0.95

##

standard_error <- population_sd / sqrt(sample_size)

confidence_onesided <- 1 - (1-confidence)/2
t_crit <- qt(confidence_onesided, df=sample_size-1)

# Create data frame for plotting
x <- seq(population_mean - 4*standard_error, population_mean + 4*standard_error, length.out = 1000)
y <- dt((x - population_mean) / standard_error, df=sample_size-1)

ci_min <- population_mean - t_crit * standard_error
ci_max <- population_mean + t_crit * standard_error

# Let's simulate an experiment where we repeat the measurement 200 times
# The confidence intervals overlap around the true mean,
# but there are some that are further away.

n_samples <- 200  # Number of samples to take

samples_df <- data.frame(
  sample_id = 1:n_samples,
  sample_mean = numeric(n_samples),
  ci_min = numeric(n_samples),
  ci_max = numeric(n_samples)
)

for(i in 1:n_samples) {
  sample_data <- rnorm(sample_size, mean=population_mean, sd=population_sd)
  samples_df$sample_mean[i] <- mean(sample_data)
  sample_se <- sd(sample_data) / sqrt(sample_size)
  samples_df$ci_min[i] <- samples_df$sample_mean[i] - t_crit * sample_se
  samples_df$ci_max[i] <- samples_df$sample_mean[i] + t_crit * sample_se
}

For a one-sample t-test, rejecting the null hypothesis with a p-value of 0.05 is equivalent of checking whether the mean of the null hypothesis is contained in the confidence interval of our data.

Sample size and power

We should be convinced that the sample size is crucial to determine how sensitive our test is to differences in the means. A large sample size will make the standard error of our estimated mean very small, regardless of the inherent variability of the population.

In many studies, we use inferential statistics to make a decision: either reject the null hypothesis, or fail to reject it.

Remember! You cannot “accept” the null hypothesis in frequentist statistics. A high p-value can either mean that:

  • the null hypothesis is true; or
  • the evidence is insufficient

The correct way of presenting a result where the p-value is >0.05 (or other specified value) is to say the evidence does not support rejection of the null hypothesis

If we want to make a binary decision (true/false), but our p-value is a continuous variable, we have to choose a significance level (usually termed \(\alpha\)), a value below which we consider it to be an “acceptable risk” to reject the null hypothesis. Whatever value we choose, we have to remember that if we repeat an experiment enough times, by definition the test will reject the null hypothesis \(\alpha\) percent of the times!

Note: The significance level is commonly set to 0.05 (5%) in many studies. There is no mathematical reason for this, but it is a design choice. You can choose whether 5% is an acceptable risk, to reject the null hypothesis when in fact it should not be rejected. If the null hypothesis is “this intervention will kill you”, you might want to be more careful when rejecting it.

Moreover, a 5% significance level means that 1 in 20 studies will demonstrate an effect when there is none.

Effect size, \(\alpha\) and \(\beta\) errors

When designing a study, we need to make some decisions. These decisions are dictated by our own knowledge of the problem, not by statistics. The first decision we need to take is what minimum effect size we want to detect. The minimum effect size is how big the difference between our studied variable (e.g. the means) and the null hypothesis needs to be to be able to be detected by my study. If we are studying a weight loss drug, for example, do we care about an average loss of 100 grams? Or 1 kg? Or 10 kg? This effect is often expressed in a normalized way, in relation to the (estimated) standard deviation of the population, because it makes effects comparable with each other. This way of expressing an effect size is called Cohen’s d:

\[d=\frac{\bar x_1 - \bar x_2}{s}\] where s is the pooled standard deviation of the two samples:

\[s=\sqrt{\frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}}\] The second decision is the acceptable risk that we reject the null hypothesis when it should not be rejected (false positive, type I error). This is the significance level \(\alpha\) that we already talked about.

The third decision is the acceptable risk that we do not reject the null hypothesis when in fact it should have been rejected (false negative, type II error). This is called \(\beta\) and it is commonly set to a value of 0.20. We are usually willing to accept more false negatives than false positives because, as explained before, the risk is that the study is inconclusive, since we cannot accept the null hypothesis.

The probability of rejecting the null hypothesis when it is in fact false (true positive result) is calculated as \(1-\beta\) and is called the statistical power of the study.

Below is an example on how these quantities relate to each other.

## Parameters

standard_deviation <- 3
effect_size <- 0.5
sample_size <- 26
alpha <- 0.05

##

mean2 <- effect_size*standard_deviation
standard_error <- standard_deviation / sqrt(sample_size)

x <- seq(-10, 10, length.out = 1000)
y1 <- dnorm(x, mean=0, sd=standard_deviation)
y2 <- dnorm(x, mean=mean2, sd=standard_deviation)

t1 <- dt(x/standard_error, df=sample_size-1)
t2 <- dt((x-mean2)/standard_error, df=sample_size-1)
significance_threshold <- qt(1-alpha, df=sample_size-1)*standard_error
# difference between the mean2 and the significance threshold for the calculation of the TypeII error
t_stat_threshold <- (significance_threshold - mean2)/standard_error
beta <- pt(t_stat_threshold, df=sample_size-1)
power <- 1 - beta

df = data.frame(
  x=x, y1=y1, y2=y2, t1=t1, t2=t2
)

densities <- stack(df[c("t1", "t2")])
df_long <- data.frame(
  x = rep(df$x, 2),
  distribution = densities$ind,
  density = densities$values
)

# Create the plot
ggplot(df_long, aes(x = x, y = density, color = distribution)) +
  geom_line() +
  geom_area(data = subset(df_long, 
                         distribution == "t1" & 
                         x >= significance_threshold),
            aes(x = x, y = density),
            fill = "red", alpha = 0.3) +
  # Fill area for t2 (red curve) below threshold
  geom_area(data = subset(df_long, 
                         distribution == "t2" & 
                         x <= significance_threshold),
            aes(x = x, y = density),
            fill = "green", alpha = 0.3) +
  scale_color_manual(values = c("t1" = "blue", "t2" = "red")) +
  theme_minimal() +
  xlim(-5,5) + 
  labs(
    title = bquote("Type I and type II errors ("*alpha*"="*.(round(alpha,3))*", "*beta*"="*.(round(beta,2))*", power="*.(round(power,2))*")"),
    subtitle = "The area shaded in red corresponds to the Type I error, the area shaded in green corresponds to the Type II error",
    color = "Distribution"
  )

Tests don’t replace your reasoning

We have now explored the t-test and in general mean testing in great detail. R has some very convenient functions to calculate t tests, without needing to manually compute the statistics.

Let’s make an example.

# Generate our mystery distributions
n <- 5000

# Magic happens here
summary_data <- data.frame(
  Distribution = c("dist_A", "dist_B"),
  Mean = c(mean(dist_A), mean(dist_B)),
  `Standard Deviation` = c(sd(dist_A), sd(dist_B))
)

knitr::kable(summary_data)
Distribution Mean Standard.Deviation
dist_A 0 2.195288
dist_B 0 2.195288 
# run a t-test. This gives us a nice summary already
t.test(dist_A, dist_B)
## 
##  Welch Two Sample t-test
## 
## data:  dist_A and dist_B
## t = -1.2904e-15, df = 9998, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.08606412  0.08606412
## sample estimates:
##     mean of x     mean of y 
## -2.832695e-17  2.832695e-17

The t-test is giving us a lot of useful information. The p-value is one, which is not surprising since both distributions have the same mean (0). But can we say that the two distributions are equal?

ggplot() +
  geom_histogram(aes(x=dist_A), bins = 200, fill="blue", alpha=0.3) +
  geom_histogram(aes(x=dist_B), bins = 200, fill="red", alpha=0.3) +
  xlim(-20,20) +
  theme_minimal()

These two distributions are clearly not equal. This is not a problem of normality of the source distributions, but rather that the mean of both distribution is effectively equal. These distributions are skewed, and the mean is not necessarily a good metric to compare them.

Important! A t-test gives us the correct answer, saying that there is no evidence to say that the distributions have different means, but this is the wrong question! Always check the plots of the distributions. Normality is not a requirement, but, in case of means, symmetry is.

A more robust approach

A potentially more robust approach is to use a test whose null hypothesis is of equal medians. In case of symmetric distributions, this hypothesis is equivalent, but in case of skewed distributions, or in case of outliers, the median is often a more useful metric.

Such test is the Mann–Whitney U test (also called the Mann–Whitney–Wilcoxon (MWW/MWU), Wilcoxon rank-sum test, or Wilcoxon–Mann–Whitney test).

In R, it is called wilcox, and we can calculate this test with the following code:

# we need to explicitly indicate that we want confidence intervals
wilcox.test(dist_A, dist_B, conf.int=TRUE, conf.level=0.95)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  dist_A and dist_B
## W = 8533583, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -0.9397100 -0.8328409
## sample estimates:
## difference in location 
##             -0.8868741

The test gives us a “location shift” parameter. This parameter is the median difference between two randomly selected observations. This assumes that the two distributions have the same shape, but shifted in medians. In case of distributions with different shapes (as the ones here, where they are mirrored), it might not be equal to the difference in medians. The confidence intervals are for this location shift parameter.

The full summary of the two distributions also tells us that, despite having equal mean, the median and the quartiles are very different:

summary(dist_A)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.6588 -1.1710 -0.6898  0.0000  0.3345 29.2743
summary(dist_B)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -29.2743  -0.3345   0.6898   0.0000   1.1710   1.6588

Multiple Testing and P-value Distribution

We mentioned earlier that, even when there is no effect, a test with a significance level of 0.05 will turn out to be true once every 20 times. Let’s demonstrate this.

## Parameters

n_tests <- 5000
n_samples <- 40
significance_level <- 0.05

##
p_values_t <- numeric(n_tests)
p_values_wilcox <- numeric(n_tests)

for(i in 1:n_tests) {
  x <- rnorm(n_samples)
  y <- rnorm(n_samples)
  p_values_t[i] <- t.test(x, y)$p.value
  p_values_wilcox[i] <- wilcox.test(x, y)$p.value
}

n_bins <- 50

# Plot p-value distribution
t_plot <- ggplot(data.frame(p_values = p_values_t), aes(x = p_values)) +
  geom_histogram(bins = n_bins, fill = "steelblue", alpha = 0.7) +
   geom_hline(yintercept = n_tests/n_bins, color="green") + 
  geom_vline(xintercept = significance_level, color = "red", linetype = "dashed") +
  labs(x = "P-value", y = "Count",
       title = "Distribution of P-values Under H0 (t-test)",
       subtitle = "Red line: significance level, Green line: uniform distribution") +
  theme_minimal()

wilcox_plot <- ggplot(data.frame(p_values = p_values_wilcox), aes(x = p_values)) +
  geom_histogram(bins = n_bins, fill = "steelblue", alpha = 0.7) +
   geom_hline(yintercept = n_tests/n_bins, color="green") + 
  geom_vline(xintercept = significance_level, color = "red", linetype = "dashed") +
  labs(x = "P-value", y = "Count",
       title = "Distribution of P-values Under H0 (ranksum)",
       subtitle = "Red line: significance level, Green line: uniform distribution") +
  theme_minimal()

grid.arrange(t_plot, wilcox_plot, ncol=2)

# Calculate proportion of false positives
sprintf("Proportion of p < 0.05 (t-test): %.3f", mean(p_values_t < significance_level))
## [1] "Proportion of p < 0.05 (t-test): 0.049"
sprintf("Proportion of p < 0.05 (ranksum): %.3f", mean(p_values_wilcox < significance_level))
## [1] "Proportion of p < 0.05 (ranksum): 0.050"

This is a very important finding to keep in mind. If a study contains multiple comparisons, it is important to correct for this statistical fluctuation.

When should we correct for multiple comparisons?

  • When we run multiple tests, and we do not care about each one specifically, but we want only to report the “significant” ones.
  • When we run the same experiment multiple times.

The most simple way of correcting for multiple comparisons is to use the Bonferroni correction. This is the most conservative approach, and consists in lowering the significance threshold by dividing it by the number of comparisons.

P-value Distribution with High Power

It is also important to observe what happens when there is a true effect:

## Parameters

effect_size <- 0.9

##

p_values_effect <- numeric(n_tests)

for(i in 1:n_tests) {
  x <- rnorm(n_samples)
  y <- rnorm(n_samples, mean = effect_size)
  p_values_effect[i] <- t.test(x, y)$p.value
}

n_bins <- 50

# Plot
ggplot(data.frame(p_values = p_values_effect), aes(x = p_values)) +
  geom_histogram(bins = n_bins, fill = "steelblue", alpha = 0.7) +
  geom_vline(xintercept = significance_level, color = "red", linetype = "dashed") +
  geom_hline(yintercept = n_tests/n_bins, color="green") + 
  labs(x = "P-value", y = "Count",
       title = paste("P-value Distribution with True Effect (Effect Size =", effect_size, ")"),
       subtitle = "The red line shows the significance level, the green line shows the uniform distribution under H0") +
  theme_minimal()

Important note! With a high-powered test (large sample size, large effect), observing p-values close to 0.05 (like 0.04) is actually more probable under the null hypothesis than under a true effect! Beware of large studies with p-values close to the significance threshold!

Key Takeaways

  1. Always be explicit in your hypotheses. Justify the choices you make for your sample sizes, effect sizes, significance thresholds.
  2. These are study design parameters. You choose them. The values you choose are not dictated by statistics, but by your knowledge of the problem.
  3. Be thorough in your reporting. Never say p < 0.05 or p > 0.05. Always report actual p-values, effect sizes, confidence intervals.
  4. Understand the assumptions of the tests. Testing for mean has a specific meaning. Nonparametric tests like the rank sum are generally preferable.
  5. A very small p-value with a large sample size doesn’t necessarily mean a clinically significant effect.
  6. A p-value close to the significance threshold with a large sample size might not be strong evidence of an effect.