What is a p-value?

In order to understand what a p-value is and what it measures, we need to go deeper into the theory behind them (the figure below illustrates what is discussed in this article). P-values belong to an statistical paradigm known as frequentist statistics (this is what you were taught in high school and basic stats courses). When analysing data using this paradigm, one needs to assume a “hypothetical infinite population, of which the actual data are regarded as constituting a random sample” (Fisher, 1923). Whether that hypothetical population does exist or not is irrelevant, as the ultimate role of this random population is to facilitate making statements about phenomena in the presence of unexplained variance or noise. Of course, this approach is most useful when there is an actual population from which the data were sampled (or when individuals were randomly assigned to groups), as it becomes easier to find a proper distribution to describe the population.

Since the data are interpreted as a random sample, any quantity calculated from the data (e.g. the arithmetic mean of observed values, \(\bar{x}\)) can be treated as a random variable. That is, each random sample from the population will yield a different value of the same quantity (e.g. different sample means). Also, it means that these quantities will follow a probability distribution. For example, if we take multiple random samples from the same population, and calculate the mean of each of these samples, some values of the sample mean will be more common that others. This distribution is what we call the sampling distribution and the quantity derived from samples is called statistic. In frequentist statistics, any probability statement (e.g. p-values, confidence intervals) always refers to a statistic and its associated sampling distribution.

The other important concept behind p-value is that of hypothesis testing. A null hypohtesis is a statement about differences between groups or phenomena. The null hypothesis represnets a default skeptical position (i.e. no difference). The goal of null hypothesis significant testing (NHST) is to challenge the null hypothesis given observed data. This procedure requires building a statistical model to describe the random population from which the data could have been sampled, assuming that the null hypothesis is true. Then, one calculates how likely it is for the observed data to occur assuming it was sampled from that population. The more unlikely the data are, the less likely the null hypothesis becomes. However, comparing the observed data to all possible alternative data is too difficult, so the problem is often simplified to comparing a statistic calculated from the data with all the possible values of the same statistic that could have been obtained from the population. In practice, this means that we need to construct a sampling distribution and evaluate how unlikely the observed statistic is within that sampling distribution.

The final complication is that the quantities that we use as statistic are often continuous measures of the data (even when the data is composed of discrete values) but the basic laws of probability tells us that the probability of any particular value of a continuous random variable is exactly zero (there are many resources online explaining this idea, for example, this one). Therefore, calculating the probability of a particular statistic value in a sampling distribution is futile. As alternative, one can calculate the probability of all values more extreme than the observed (i.e. the shaded area in the figure above). The smaller this probability, the more unlikely the value is. And this probability is the p-value!

Since all this was very abstract, we will go through a classic example: the one sample t-test. I will illustrate with a computational approach how the sampling distribution arises from the population and how it relates to the sample statistic. I chose R as the programming language for this article as this is a common language for statistics, but the same approach could have done in Matlab/Octave, Python or Julia.

Example for the one sample mean test (one sample t-test)

Let’s imagine that we have a small dataset such as:

sample = c(1.52, 5.24, -0.23, 2.47, 2.63)

Most values are positive and we may want to test whether the underlying effect is indeed positive. We will assume that this data could be treated as a random sample from a Normal population with a positive mean (i.e. \(\mu > 0\)). As described above, we set as null hypothesis (\(H_0\)) the more conservative option, which in the case it would be \(\mu = 0\), while the contesting or alternative hypothesis (\(H_1\)) represents positive case (\(\mu > 0\)). Formally, this is written as:

\[ \begin{align} H_{0} &: \mu = \mu_0 \\ H_{1} &: \mu > \mu_0 \end{align} \]

Where \(\mu_0 = 0\) in this case. Assuming that the sampling procedure ensures statistical independence and given the assumptions in the above, we could use the t statistic, \((\bar{x} - \mu_0)/(s/\sqrt n )\), where \(\bar{x}\) is the sample mean, \(s\) is the standard deviation of the sample and \(n\) is the sample size. The reason for using this statistic is that the sampling distribution is known to be a Student’s t distribution. The proof for this was given by William Gosset in 1908 (signed with the pseudonym Student). If you want to check what it takes to make such a proof, take a look at the original paper.

Rather than going through the mathematical proof, I will use a simulation-based procedure that I find more intuitive and easier to follow, and it is also very general, while producing the same results. Effectively, we will construct a Monte Carlo approximation to the sampling distribution of the t statistic. The computational procedure is as follows:

1. Estimate all the parameters of the population from the sample, except for the parameter that is already fixed by the null hypothesis. In this example, this means that we estimate the standard deviation of the Normal population (\(\sigma\)) from the sample, while fixing the mean (\(\mu\)) to 0.

2. Generate \(N\) random samples from the population of size \(n\) using the Monte Carlo method. In this case \(n\) is equal to 5 and \(N\) can be as large as you want (the larger the better).

3. For each sample generated in step 2, calculate the statistic. In this this case, we calculate \((\bar{x} - \mu_0)/(s/\sqrt n )\) for each sample.

4. The \(N\) statistics calculated in step 3 is a large sample of values from the sampling distribution. The p-value can then be calculated as the fraction of this sample that exceeds (or is lower, depending on the tail being tested) than the value of the statistic in the actual data.

for(lib in c("furrr", "distr6")) {
  library(lib, character.only = TRUE)
}

plan(multiprocess)

set.seed(2019)

sigma_hat = sd(sample)

mu = 0

population = Normal$new(mean = mu, sd = sigma_hat)

calc_statistic = function(x, mu0 = 0) {
  (mean(x) - mu0)/(sd(x)/sqrt(length(x)))
}

N = 5e4
sampling_distribution = future_map_dbl(1:N, ~ calc_statistic(population$rand(5)))

plot(density(sampling_distribution,), xlab = "t statistic", main = "", xlim = c(-10,10)) # Histogram of sampling_distribution
curve(dt(x,4), -10, 10, add = T, col = 2, n = 1e3) # Theoretical distribution t(n-1)
abline(v = calc_statistic(sample), col = 3) # Observed statistic

sum(sampling_distribution >= calc_statistic(sample))/N

## [1] 0.02874

1 - pt(calc_statistic(sample), 4)

## [1] 0.02946129

t.test(sample, alternative = "greater")$p.val

## [1] 0.02946129

N = 5e4
sampling_distribution = future_map_dbl(1:N, ~ mean(population$rand(5)))

plot(density(sampling_distribution, bw = "SJ"), main = "", xlab = "sample mean") # Histogram of sampling_distribution
abline(v = mean(sample), col = 2) # Observed statistic

sum(sampling_distribution >= mean(sample))/N

## [1] 0.0043

Final remarks

As has been shown in this article, in order to calculate the p-value of an statistic we need to know its sampling distribution, which depends on the model assumed for the population. Traditionally, one had to derive the sampling distribution mathematically (if possible) for every combination of statistic and population model. This greatly restricts the number of scenarios which can be handled by traditional statistical tests, which leads to that dreaded question: “Which statistical test should I use for my data?”. There are dozens of decision trees that the practioner of statistics is supposed to follow (you can see some in this paper).

Nowadays, you can use a simulation approach such as the one presented in this article. With such an approach, you can always calculate the p-value for any statistic and population model, even if no one has ever tried that particular combination before. Restricting oneself to analytical methods made sense in the early 20th century when many of these tests were developed, but in the 21st century you can use computers! In the next post I will show how to apply this method to real world data that was not (and could not be) normally distributed. Stay tuned!