P-value is the probability that the treatment effect is larger than zero (under certain conditions)

A quick t-test refresher

If we were to evaluate our results using a t-test, we would:

• Compute sample means: $\bar{Y_A} = \text{mean}(Y_A)$, $\bar{Y_B} = \text{mean}(Y_B)$

• Compute sample variances: $s_A^2\ =\text{var}(Y_A)$, $s_B^2\ =\text{var}(Y_B)$

Then, enabled by CLT, we can claim that the difference of the means between the two groups is normally distributed (that distribution is referred to as the “sampling distribution”), with the following parameters:

• Mean: $\bar{Z} = \bar{Y_A} - \bar{Y_B}$

• Standard error: $s_Z = \sqrt{\frac{s_A^2}{n_A} + \frac{s_B^2}{n_B}}$

Finally, we calculate the test statistic $t = \bar{Z} / s_Z$, and look up the associated CDF of the sampling distribution at that point (a.k.a p-value).

This is how it may look numerically.

set.seed(42)

n = 6000

effect_size = 0.05



A = rnorm(n, 0, 1)

B = rnorm(n, effect_size, 1)



sampling_distribution = function(X, Y) {

  list(

    mean = mean(X) - mean(Y),

    sd = sqrt(var(X) / n + var(Y) / n)

  )

}



sdist = sampling_distribution(B, A)

t = sdist$mean / sdist$sd



as_tibble(sdist) |> 

  pivot_longer(everything(), names_to = "statistic") |> 

  kbl(caption='Summary statistics of the sampling distribution') |> 

  kable_minimal(full_width = T, position='left')

tibble(

  statistic = c('Test statistic', 'P-value'),

  value = c(t, (1 - pnorm(t)) * 2)

) |> kbl(caption='Hypothesis test outcomes') |> 

  kable_minimal(full_width = T, position='left')

As a visual reminder - the p-value is twice the CDF of the sampling distribution in the range $[-\infty, 0]$ (twice, because we’re testing for a two-sided alternative hypothesis).

Mathematically, we don’t need to calculate the test statistic. We can get the CDF directly using the parameters of our sampling distribution.

#this is also a p-value!

format(pnorm(0, mean = sdist$mean, sd = sdist$sd, lower.tail=T) * 2, digits=4)

[1] "0.0224"

We can double-check using the built-in R functions.

t.test(B, A)

    Welch Two Sample t-test



data:  B and A

t = 2.2835, df = 11997, p-value = 0.02242

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

 0.00593807 0.07794067

sample estimates:

   mean of x    mean of y 

 0.033779917 -0.008159451 

Great, we can conclude that, if our hypothesis of treatment having no effect were true, we would be observing a difference between two groups at least as large as we did (we observed a difference of 0.041) roughly 2% of the time.

Sounds rare enough that we could reject our null hypothesis and conclude that there indeed is a difference.

We also conclude that the true effect should be somewhere between 0.005 and 0.078 95% of the time if we were to repeat this experiment an infinite number of times. As for this particular run - well, it’s either inside that range or not. T-tests (or frequentist approaches, more generally) don’t tell you anything about this particular experiment. Just long-term probabilities.

A mouthful. But hey, this isn’t a post about issues with null hypothesis testing. I wrote that one already.

A Bayesian approach

Generally, saying “I’m using a Bayesian approach” is like saying “I’m using linear algebra”. It’s not a statistical method. We need an estimand! We will use the simplest one: the probability that the treatment is larger than 0.

Next, we will need uninformative priors. We can achieve that by choosing normal distribution priors with variance approaching infinity (making any number in the real domain roughly equally likely), i.e.:$Y_A, Y_B \sim N(0, 1/\tau)$, where $\tau$ is some tiny number, e.g. 0.001.

Note that we did not choose a prior for the treatment effect itself. Instead, we’re effectively saying that both experiment arms have the same mean and can take almost any value with the same likelihood.

However, we can see what kind of treatment effect prior distribution our choice implies by leveraging the fact that a difference between two independent normal distributions is just another normal distribution:

$$\text{if } X \sim N(\mu_x, s_x^2), Y \sim N(\mu_y, s_y^2)\, \text{ then } (X - Y) \sim N(\mu_x - \mu_y, s_x^2 + s_y^2)$$

difference_of_two_normals = function(X, Y) {

  list (

    mean = X$mean - Y$mean,

    tau = 1 / ((1 / X$tau) + (1 / Y$tau)),

    sd = sqrt((1 / X$tau) + (1 / Y$tau))

  )

}





difference_of_two_normals(list(mean=0, tau=0.001), list(mean=0, tau=0.001)) |> 

  as_tibble() |> pivot_longer(everything(), names_to = 'statistic') |>

  kbl(caption='Implied treatment effect prior distribution') |> 

  kable_minimal(full_width = T, position='left')

We can see that choosing a prior for experiment results of zero mean and $\tau=0.001$ (or, equivalently, the standard deviation of $\sqrt(\frac{1}{0.001})=31.6$) yields a prior on the treatment effect that also has a zero mean but has even more uncertainty, with a standard deviation of $44.7$.

Next, we need to get to posterior distributions. Luckily, a normal prior and a normal likelihood form a conjugate pair that yields a normal posterior distribution, so we can do this algebraically.

Specifically, suppose the prior distribution is $P_0 \sim N(\mu_0, \tau_0^{-1})$, the observed mean is $\bar{X}$ and observed precision is $\bar{\tau} = \frac{1}{s^2_x}$. Then the posterior is $P \sim N \left(\frac{\tau_0 \mu_0+\bar{\tau} \bar{X} n}{\tau_0+ n \bar{\tau}},\left(\tau_0+n \bar{\tau}\right)^{-1}\right)$. Let’s calculate the posteriors for our samples.

calculate_normal_posterior = function(X, prior) {

  n = length(X)

  observed_mu = mean(X)

  observed_tau = 1 / var(X)

  

  list(

    mean = (prior$tau * prior$mean + observed_tau * observed_mu * n) / (prior$tau + n * observed_tau),

    tau = (prior$tau + n * observed_tau),

    sd = sqrt(1 / (prior$tau + n * observed_tau))

  )

}



A_posterior = calculate_normal_posterior(A, list(mean=0, tau=0.001))

B_posterior = calculate_normal_posterior(B, list(mean=0, tau=0.001))



bind_rows(A_posterior, B_posterior, .id = 'ID') |> 

  mutate(ID = str_replace_all(ID, c('1'='A', '2'='B'))) |>

  kbl(caption='Posterior distributions of experiment arms') |> 

  kable_minimal(full_width = T, position='left')

We have posterior distributions for the data of the two samples, but we are after the treatment effect - i.e., the difference between the two samples. We can use the same approach as we did previously when estimating the implied treatment effect prior previously and calculate it, as the two posteriors can be assumed to be independent, too.

treatment_posterior = difference_of_two_normals(B_posterior, A_posterior)

treatment_posterior |> as_tibble() |> pivot_longer(everything(), names_to = 'statistic') |>

  kbl(caption='Posterior distribution of the treatment effect') |> 

  kable_minimal(full_width = T, position='left')

Let’s visualize it.

ggplot() + geom_function(

  fun = dnorm, 

  args = within(treatment_posterior, rm("tau"))

) + xlim(-0.1, 0.1) + 

  geom_vline(xintercept=0, color='red', linetype='dashed') + 

  xlab("") + ggtitle("PDF the treatment effect") + ylab("Density") +

  theme_light() + theme(panel.grid=element_blank())

That looks familiar! Let’s overlay the sampling distribution from the t-test:

ggplot() + geom_function(

  fun = dnorm, 

  args = within(treatment_posterior, rm("tau")),

  aes(color = 'Posterior of a treatment effect'),

  linetype='dashed'

) + xlim(-0.1, 0.1) + 

  geom_function(

  fun = dnorm, 

  args = sdist,

  aes(color = 'Sampling distribution'),

  linetype='twodash'

) + xlab("") + ggtitle("Treatment effect posterior distribution vs. Sampling distribution") + 

  ylab("Density") + labs(color='') +

  theme_light() + theme(panel.grid=element_blank(), legend.position='bottom')

They are identical! This means that the probability of the treatment effect being smaller than zero is—you guessed it—equal to 1/2 of a two-sided p-value—or, in other words, identical to the one-sided p-value.