Aurimas Online https://aurimas.eu/blog Sun, 09 Jun 2024 20:12:32 +0000 en-US hourly 1 https://wordpress.org/?v=7.0 161723146 P-value is the probability that the treatment effect is larger than zero (under certain conditions) https://aurimas.eu/blog/2024/06/p-value-is-the-probability-that-the-treatment-effect-is-larger-than-zero-under-certain-conditions/ Sun, 09 Jun 2024 20:12:32 +0000 https://aurimas.eu/blog/?p=635 All code behind this blog post is available on GitHub.

P-value is the probability that treatment effect is larger than zero. Under certain conditions.

Ah, the beautiful concept of p-value. The probability of observing the result as extreme or more extreme than the data at hand, assuming the null hypothesis is true. Not a single person ever simplified it to something else. Like.. I don’t know.. the probability that the null hypothesis is true? Or, in an A/B test, the probability that treatment is greater/smaller than zero? 👀

Doing that is equivalent to committing an original data science sin. You’d be banished to the ungodly world full of local optima, spaghetti SQL, Excel downloads and out-of-date dashboards. Or simply wouldn’t pass an screening round in a hiring process.

BUT.

P-value can be equal to the probability that the treatment effect is greater/smaller than zero. In the context of A/B testing, where we are in a straightforward univariate hypothesis testing setting, it is as simple as employing a Bayesian approach with an uninformative prior.

This post, I admit, is partly for kicks. But it has some practical implications, too:

  • Most articles on using Bayesian A/B testing will tell you that using uninformative priors is not a very good idea. While I agree, the fact that using them yields results equivalent to a most popular frequentist procedure means that they aren’t a terrible choice. There are some vendors out there that scaremonger into shying away from using Bayesian approaches. This post shows that if you choose to go with uninformative priors, you’ll be shooting yourself into the foot as much as you would with a simple t-test.

  • The bad thing, however, is that p-values have such an interpretation under very unrealistic assumptions of prior distributions. In practice, thus, interpreting p-value as a direct measure of evidence or using uninformative priors is a bad idea.

The Gaussian case

Let’s consider the following setup:

  • We’re interested in measuring the impact of a simple A/B test.

  • Our outcome metric happens to be normally distributed: YA∼N(μA,σA)Y_A \sim N(\mu_A, \sigma_A) and YB∼N(μB,σB)Y_B \sim N(\mu_B, \sigma_B). For simplicity, we will set μA=0,μB=0.05, and σA,σB=1\mu_A = 0, \mu_B = 0.05, \text{ and } \sigma_A, \sigma_B = 1. In other words, equal variance, and a treatment effect of 0.050.05.

Suppose we ran an experiment and collected 6,000 samples in each experiment arm (that’s 78% power, assuming we set MDE equal to the true unobservable population effect).

A quick t-test refresher

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

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

  • Compute sample variances: sA2=var(YA)s_A^2\ =\text{var}(Y_A), sB2=var(YB)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: Z‾=YA‾−YB‾\bar{Z} = \bar{Y_A} - \bar{Y_B}

  • Standard error: sZ=sA2nA+sB2nBs_Z = \sqrt{\frac{s_A^2}{n_A} + \frac{s_B^2}{n_B}}

Finally, we calculate the test statistic t=Z‾/sZt = \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.

Show code 
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')
Summary statistics of the sampling distribution
statistic value
mean 0.0419394
sd 0.0183665
Show code 
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')
Hypothesis test outcomes
statistic value
Test statistic 2.2834717
P-value 0.0224026

As a visual reminder - the p-value is twice the CDF of the sampling distribution in the range [−∞,0][-\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.

Show code 
#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.

Show code 
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.:YA,YB∼N(0,1/τ)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:

if X∼N(μx,sx2),Y∼N(μy,sy2) then (X−Y)∼N(μx−μy,sx2+sy2) \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)

Show code 
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')
Implied treatment effect prior distribution
statistic value
mean 0.00000
tau 0.00050
sd 44.72136

We can see that choosing a prior for experiment results of zero mean and τ=0.001\tau=0.001 (or, equivalently, the standard deviation of (10.001)=31.6\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.744.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 P0∼N(μ0,τ0−1)P_0 \sim N(\mu_0, \tau_0^{-1}), the observed mean is X‾\bar{X} and observed precision is τ‾=1sx2\bar{\tau} = \frac{1}{s^2_x}. Then the posterior is P∼N(τ0μ0+τ‾X‾nτ0+nτ‾,(τ0+nτ‾)−1)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.

Show code 
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')
Posterior distributions of experiment arms
ID mean tau sd
A -0.0081594 5873.761 0.0130479
B 0.0337799 5985.180 0.0129259

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.

Show code 
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')
Posterior distribution of the treatment effect
statistic value
mean 0.0419394
tau 2964.4734513
sd 0.0183665

Let’s visualize it.

Show code 
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:

Show code 
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.

This is no coincidence. We can prove that mathematically. Here’s the posterior formula again:

P∼N(τ0μ0+τ‾X‾nτ0+nτ‾,(τ0+nτ‾)−1) 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)

Now, if we let the prior variance go to infinity, the τ0\tau_0 term goes to zero. In that case, the posterior simplifies to:

P∼N(X‾,(nτ‾)−1) P \sim N \left(\bar{X},\left(n \bar\tau\right)^{-1}\right)

We can plug it into the formula for calculating the difference between two normal distributions:

(X−Y)∼N(μx−μy,sx2+sy2)=N(X‾−Y‾,1nXτX+1nYτY)=N(X‾−Y‾,sX2nX+sY2nY) (X - Y) \sim N(\mu_x - \mu_y, s_x^2 + s_y^2) = N(\bar{X} - \bar{Y}, \frac{1}{n_X \tau_X} + \frac{1}{n_Y \tau_Y}) = N(\bar{X} - \bar{Y}, \frac{s_X^2}{n_X} + \frac{s^2_Y}{n_Y})

What we get is precisely the sampling distribution used in the frequentist t-test.

So, no, it’s not true that using uninformative priors will yield rubbish results in an A/B test. In fact, you’ll get precisely the same results as a one-sided p-value, which also means you have the same statistical guarantees. If your decision criterion is to look at experiment results once and adopt the treatment when the probability that it is above zero is 95% or more (and thus the probability it’s below zero is 5% or less), you will be right 95% of the time, just like you would if you were to rely on p-values and statistical significance.

It also shows that Bayesian methods are not immune to peeking (as the results are equivalent to a t-test which isn’t either). Peeking is either addressed by informative priors or by reframing the errors you care about altogether.

One-sided vs. two-sided p-values

At first, the fact that the Bayesian approach leads to a one-sided p-value confused me a bit. I always found one-sided tests weird. What do you mean you are testing only in one direction? I understand that there may be situations where only a single direction matters, but in an A/B testing context, where the outcome metric is a business KPI, surely you want to know if it’s positive or negative? And how stupid it would be when you upfront choose to test only in a positive direction but the observed effect is negative and then you don’t report whether it’s statistically significant from zero?

For a minute, I thought perhaps all those people arguing for one-sided tests were right. But then I came across this thread on X/Twitter:

And indeed, technically, if you do a two-sided t-test, all you obtain is an estimation of whether the effect is not zero. To be able to claim that it is larger/smaller than zero, you need two one-sided t-tests instead. Of course, because you are now doing two tests instead of one, you need to adjust for multiple comparisons, so to achieve a 5% Type I error rate, you need to use α=0.025\alpha=0.025.

This is all horribly confusing if you ask me. Let’s move on.

A non-Gaussian case

You may rightfully wonder if the above result extends to non-Gaussian situations. I can’t prove it mathematically, but numerical simulations tell me the answer is yes. Let’s consider the most popular A/B testing situation involving binary outcomes (e.g., conversion rates).

  • Our outcome metric follows the Bernoulli distribution: X∼Bernoulli(p)X \sim \text{Bernoulli}(p). Because the individual user observations are independent, we can interpret the entire sample with nn observations and kk successes as a binomial distribution with rate pp.

  • We will use uninformative beta distribution priors P0∼Beta(α0=1,β0=1)P_0 \sim \text{Beta}(\alpha_0=1,\beta_0=1)

  • The posterior of a beta prior and binomial likelihood is another conjugate pair, resulting in a beta distribution with parameters α=α0+k,β=β0+(n−k)\alpha=\alpha_0 + k, \beta = \beta_0 + (n - k).

We will also need to be able to calculate a difference between two random variables following Beta distributions. A closed-form solution exists to calculate this difference, but it involves Appell F1 hypergeometric functions that I couldn’t get to work with large nn. So instead, we will rely on the fact that a sum/difference of two independent PDFs can be calculated using convolutions.

Let’s generate some data. We’ll use a baseline of 0.170.17, an effect size of 0.020.02, and a group size of 6,000, which should yield us ~80% power.

Show code 
eff = 0.02

baseline = 0.17

n_group = 6000



power.prop.test(n=n_group, p1 = baseline, p2 = baseline + eff, sig.level=0.05)


     Two-sample comparison of proportions power calculation 



              n = 6000

             p1 = 0.17

             p2 = 0.19

      sig.level = 0.05

          power = 0.8137145

    alternative = two.sided



NOTE: n is number in *each* group

Let’s define a couple of helper functions.

Show code 
get_samples = function(n_group, baseline, eff) {

  list(

    A = rbinom(n_group, 1, baseline),

    B = rbinom(n_group, 1, baseline + eff)

  )

}



calculate_beta_posterior = function(X, prior){

  list(

    alpha = prior$alpha + sum(X),

    beta = prior$beta + (length(X) - sum(X))

  )

}



#calculate density at a given point X for the difference between two beta distributions

pdf_beta_diff = function(X, d1, d2) {

  

  sapply(X, function(t) {

    integrate(function(i) {

      (dbeta(i, d1$alpha, d1$beta) * dbeta(i-t, d2$alpha, d2$beta) + dbeta(i, d2$alpha, d2$beta) * dbeta(i+t, d1$alpha, d1$beta)) / 2

    }, lower=-Inf, upper=Inf, stop.on.error = F)$value

    })

  

}



#calculate a total CDF between X and Y for two beta distributions

cdf_beta_diff = function(X, Y, d1, d2) {

  integrate(pdf_beta_diff, d1, d2, lower=X, upper=Y, stop.on.error = F)$value

}

Now, let’s draw one sample and do the calculations.

Show code 
set.seed(42)

data = get_samples(n_group, baseline, eff)

sdist = sampling_distribution(data$B, data$A)

sdist |>as_tibble() |>

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

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

  kable_minimal(full_width = T, position='left')
Summary statistics of the sampling distribution
statistic value
mean 0.016500
sd 0.007012
Show code 
A_posterior = calculate_beta_posterior(data$A, list(alpha=1, beta=1))

B_posterior = calculate_beta_posterior(data$B, list(alpha=1, beta=1))

Let’s visualize the posterior vs. the sampling distribution.

Show code 
ggplot() + geom_function(

  fun = pdf_beta_diff, 

  args = list(d1=B_posterior, d2=A_posterior),

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

  linetype='twodash'

) + xlim(-0.05, 0.05) + 

  geom_function(

  fun = dnorm, 

  args = sdist,

  aes(color = 'Sampling distribution'),

  linetype='dotted'

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

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

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

And compare one-sided p-value vs. posterior probability that treatment is below zero.

Show code 
tibble(

  statistic = c('One-sided p-value', 'Probabability treatment is < 0'), 

  value = c(

    t.test(data$B, data$A, alternative='g')$p.value, 

    cdf_beta_diff(-Inf, 0, B_posterior, A_posterior)

  )

) |> kbl(caption='P-value vs. Bayesian estimate of treatment effect < 0') |> 

  kable_minimal(full_width = T, position='left')
P-value vs. Bayesian estimate of treatment effect
statistic value
One-sided p-value 0.0093164
Probabability treatment is < 0 0.0093215

Everything looks identical. To be sure, we can repeat the procedure 1000 times.

Show code 
sim_results = sapply(1:1000, function(i) {

  data = get_samples(n_group, baseline, eff)

  A_posterior = calculate_beta_posterior(data$A, list(alpha=1, beta=1))

  B_posterior = calculate_beta_posterior(data$B, list(alpha=1, beta=1))

  posterior_prob_above_zero = cdf_beta_diff(-Inf, 0, B_posterior, A_posterior)

  p_value = t.test(data$B, data$A, alternative='g')$p.value

  c(posterior_prob_above_zero, p_value)

})

I hope that’s convincing.

All those difficulties with interpreting p-values? A p-value is not a measure for evidence, etc.? Well, I think I have proven to myself (and hopefully, you, too) that if we are willing to assume that we are equally likely to observe data of any value, then under our very simple and restrictive setting of two independent samples with data generating process following Normal distribution (or distribution that converges to Normal, such as Beta):

  • A one-sided p-value is equal to the probability that the treatment is negative;

  • (1 - p-value) is equal to the probability that the treatment is positive (gasp!)

In other words, we can all run t-tests and report p-values as probabilities of treatment effects. No more worrying about explaining statistical significance!

Sure, I am joking. Or am I?

One good reason not to do that is because the assumption of flat priors is a pretty bad one. If you’re running conversion tests on a website, surely you know better than to say, “Ah, it can be anything.” So maybe that’s the answer to the people who want to misinterpret p-values as measures of evidence — “Yes, it’s a probability of treatment being better/worse, but only if you’re assuming our conversion rates can be anything.”

]]> 635 Estimating long-term detection, win, and error rates in A/B testing https://aurimas.eu/blog/2024/05/estimating-long-term-detection-win-and-error-rates-in-a-b-testing/ Mon, 27 May 2024 21:04:22 +0000 https://aurimas.eu/blog/?p=615 This blog post is also available in Quarto notebook form on GitHub.

Estimating long-term detection, win and error rates in A/B testing

Why a team that has good (but not great) ideas 75% of the time may see a win rate of just 25%, and why focusing on MDE/statistical power alone may not always be the best

Classical statistical guarantees in A/B testing (Type I and II error rates) are great at quantifying the reproducibility of a single experiment. As a reminder:

  • Type I error rate, α\alpha, guarantees that when the null hypothesis is true (typically: no difference between variants), you will only observe statistically significant results (i.e., observed data says there’s a difference) α\alpha fraction of the time. It’s the false positive rate for an experiment with no effect.

  • Type II error rate, 1−β1-\beta, guarantees that when the null hypothesis is false, you will observe data that will yield statistically significant results β\beta fraction of the time (a.k.a statistical power), conditional on sample size, effect size of interest, and chosen α\alpha rate. It’s the false negative rate for an experiment with an underlying effect.

However, they tell us little about the long-term dynamics of an “experimentation program” ( a series of experiments, each estimating the impact of a different treatment). You can’t immediately answer questions such as “What % of experiments can I expect to detect a statistically significant result in?”.

That, however, may be an interesting question to answer. Suppose you’re Head of Data Something in a company that ran 100 experiments last year, of which only 5% yielded positive statistically significant results (a bit below an oft-quoted figure in the industry). You may be wondering if that’s normal. Is it really so that 95% of product changes had zero impact? Perhaps they weren’t exactly zero, but simply smaller than you could declare as statistically significant with the sample sizes at hand?

If you work in a large company, you may find yourself at the opposite of “statistically but not practically significant.” After all, a conversion rate improvement of 0.5 percentage point may still mean many $$ annually. However, achieving 80% power on such an effect size requires a 3-4x sample size compared to an effect size of 1 percentage point. Are you willing to make that trade-off and make all experiments 3-4x slower?

Understanding the relationship between experiment-level design parameters (minimum detectable effect (MDE), sample size, power) and overall experimentation program win/detection/error rates can be pretty illuminating. So, let’s figure out how to do it.

This blog post is accompanied by an interactive Streamlit app where you can explore how all the different parameters come to life yourself.

From statistical power to detection/win rates

You may wonder why the answer to long-term detection rates isn’t simply the Type II error rate or “We will detect (i.e., declare them as statistically significant) differences when they exist 80% of the time”.

That’s because the Type II error rate is conditional on the MDE. Assuming you use the same MDE all the time, you can claim that when the actual effect is at least as large as the MDE, we will detect it at least 80% of the time.

But what if it’s not? And how exactly does one think about the “at least” part of the statement? What if actual effects are larger/smaller than MDE? Your power will be exactly 80% only if you always choose MDE precisely equal to the actual unobservable difference for that specific experiment. If you can do that, you don’t need to experiment!

Frequentist approaches tend to develop tools that find upper/lower bounds of a measure of interest, and statistical power is an excellent example of that (I saw this clearly articulated in an essay by Jacob Steinhardt titled “Beyond Bayesians and Frequentists” - I highly recommend it).

However, we are after an expectation—we want to know the average expected detection rate. For that, we’ll borrow an idea from Bayesian statistics that tends to be better suited for such purposes. Instead of working with fixed values, we will assume a distribution of likely effect sizes (a.k.a., a prior distribution) and then integrate over that distribution to get expected rates.

To be precise, by integrating over a distribution of hypothesized treatment effect sizes AA (we’ll call the distribution FAF_A), we can calculate:

  • Expected “detection rate” - the % of experiments we would declare statistically significant given they have a non-zero actual treatment effect.

  • Expected “win rate” - the % of experiments we would declare statistically significant, but we count only experiments with a positive treatment effect.

Mathematically:

expected detection rate=∫−∞+∞FA(x)Pr(reject H0|x≠0,α,n,σ2)dx \text{expected detection rate} = \int_{-\infty}^{+\infty} F_A(x) \text{Pr}(\text{reject } H_0 | x \neq 0, \alpha, n, \sigma^2) dx expected win rate=∫0+∞FA(x)Pr(reject H0|x≠0,n,α,σ2)dx \text{expected win rate} = \int_0^{+\infty} F_A(x) \text{Pr}(\text{reject } H_0 | x \neq 0, n, \alpha, \sigma^2) dx

where nn is the sample size, σ2\sigma^2 is the variance of the outcome metric, and α\alpha is the chosen statistical significance level.

One immediate thing to note is that detection/win rates do not depend on MDEMDE or β\beta (statistical power). The only levers you have to increase them are larger samples, variance reduction techniques (e.g., by using covariates, such as CUPED-like methods), better treatments (changes in FAF_A), and different α\alpha levels.

If your immediate reaction is, “OMG, a subjective assumption on effect size distribution (gasp!)”, I would like to remind you that choosing the MDE threshold for power calculations is subjective, too. Ideally, one would pick an MDE upfront and tune the sample size to achieve the desired statistical power. However, in practice, there’s always an incentive to cherry-pick an MDE itself, given a fixed sample size. It looks better to report a power of 80% with an MDE of 0.015 than a power of 50% with an MDE of 0.01.

Having said that, if you use MDE correctly and report statistical power together with the MDE, bad MDE choices are more transparent than flawed effect size distribution assumptions. But if we want to estimate average detection/win rates, that’s the cost we need to pay - there’s no free lunch!

Selecting an effect size distribution

So, how do we select a potential effect size distribution?

  • Historical experiment results could be a good starting point for understanding where most of the distribution density should be. If most experiments see effect sizes of a few percentage points, then the distribution should reflect that, too. Importantly, we should consider observed treatment effects in all experiments, not just the statistically significant ones.

  • Teams have domain expertise and are paid to do their work, so we’d hope they do better than a random flip of a coin - the mean/median of the distribution should be positive.

  • The distribution should not be symmetric. Teams may occasionally experiment with more drastic changes, but safeguards (dogfooding, small-scale pilots, UX research) should prevent ideas with a major negative impact from reaching the experimentation phase.

From a modeling perspective, a (shifted/scaled) Gamma distribution is a good option for modeling such skewed distributions with limited downside and a longer tail of positive effects.

For the sake of this blog post, I will use two distributions representing different hypothetical teams:

  • Mature product team. This team manages a mature product with a well-established user base. Experiments usually represent incremental changes and minor features, and thus, an average experiment achieves just a 0.50.5 percentage point improvement. They extensively rely on dogfooding, UXR, and other means to test ideas before they reach the experimentation phase, thus protecting the downside very well. We’ll model it with a Gamma distribution.

  • Startup product team. This team is still in the product-market fit stage. They tend to make wide-ranging changes and like moving fast, without too much research before the experimentation phase. Because it’s a less polished product, experiments tend to produce higher gains (2.02.0 percentage point improvements, on average), but they also see negative impacts from changes tested more often. We’ll model it with a Normal distribution.

Show code 
eff_size_distrs = list(

  list(

    name = 'Startup team',

    distr = 'x ~ N(0.02, 0.03)',

    generator = function() rnorm(1, mean=0.02, sd=0.03),

    density = function(x) dnorm(x, mean = 0.02, sd = 0.03) 

  ),

  list(

    name = 'Mature team',

    distr = 'x ~ Gamma(2,200)) - 0.005',

    generator = function() rgamma(1, shape=2, rate = 200) - 0.005,

    density = function(x) dgamma(x + 0.005, shape=2, rate = 200)

  )

)

Here’s how that looks visually, with some summary statistics below.

Show code 
emp = lapply(

  eff_size_distrs,

  function(d) {

    r = list()

    r[[d$name]] = sapply(1:50000, function(i) d$generator())

    r

  }

) |> bind_cols() |> pivot_longer(everything(), names_to = "Assumed distribution")



ggplot(emp) + 

  #geom_density(aes(x=value, color=`Assumed distribution`)) +

  sapply(

  eff_size_distrs, 

  function(d) stat_function(

    aes(color=d$name, fill=d$name), 

    fun = d$density,

    data = tibble(name = d$name, full_name = paste(d$name, '', d$distr)),

    geom="area",

    alpha=0.5

  )

) + 

  xlim(-0.05, 0.1) + labs(

    y='Density',

    color='Assumed effect size distribution',

    x='Effect size',

    title='Some possible effect size assumptions'

  ) + 

  theme_light() + 

  theme(

    panel.grid.major = element_blank(),

    legend.position = "none"

  ) + facet_wrap(~full_name) +

  geom_vline(xintercept = 0, linetype = 'dotted')

Show code 
emp |> group_by(`Assumed distribution`) |>

  summarize(

    `Avg. effect` = mean(value),

    `10th pctile` = quantile(value, 0.1),

    `Q50 (median)` = quantile(value, 0.5),

    `90th pctile` = quantile(value, 0.9),

    `Share of experiments with positive effect` = mean(if_else(value > 0,1,0))

  ) |> 

  arrange(-`Avg. effect`) |>

  kbl(

    caption = 'Summary statistics of assumed effect size distributions',

    digits=4

  ) |>

  kable_styling()
Summary statistics of assumed effect size distributions
Assumed distribution Avg. effect 10th pctile Q50 (median) 90th pctile Share of experiments with positive effect
Startup team 0.0201 -0.0180 0.0199 0.0584 0.7505
Mature team 0.0050 -0.0023 0.0034 0.0143 0.7366

Note that both teams are expected to produce nearly the same (~75%) share of experiments that have a positive impact. We could interpret that as both teams having a similar “quality of ideas.” It’s just that the effect sizes differ due to the different contexts in which they operate - a median experiment of the mature team is assumed to produce just a 0.3ppt effect, while the startup team - 2ppt, a 6x difference.

Estimating average expected detection/win rates

To estimate detection/win rates, we need a few more inputs. Let’s assume that:

  • Both teams are optimizing a metric that has a baseline of 0.170.17

  • Both teams run experiments with 10,000 observations per group (let’s say they both used 0.0150.015 as their MDE and chose β=0.2\beta=0.2; if you were to plug these inputs into a sample size calculation, you’d get close to 10,000).

Show code 
baseline = 0.17

required_sample_size = 10000

mde = 0.015



#check that power is ~80%

power.prop.test(n=required_sample_size, p1=baseline, p2=baseline + mde)


     Two-sample comparison of proportions power calculation 



              n = 10000

             p1 = 0.17

             p2 = 0.185

      sig.level = 0.05

          power = 0.7927863

    alternative = two.sided



NOTE: n is number in *each* group

What long-term detection rates should they expect to see?

Show code 
estimate_detection_rate = function(baseline, required_sample_size, density_func) {

  integrate(function(x) {

    p = power.prop.test(

      n=required_sample_size, 

      p1=baseline,

      p2=baseline + x

    )$power    

    density_func(x) * replace_na(p, 0)    

  }, -Inf, Inf)$value  

}



tibble(

  `Assumed effect size distribution` = sapply(eff_size_distrs, function(x) x$name),

  `P(detect | effect size distribution, N, baseline)` = sapply(

    eff_size_distrs,

    function(x) estimate_detection_rate(baseline, required_sample_size, x$density)

  )

) |> kbl(caption='Long-run detection rates', digits=2) |> kable_styling()
Long-run detection rates
Assumed effect size distribution P(detect | effect size distribution, N, baseline)
Startup team 0.78
Mature team 0.25

A massive difference! The startup team will be detecting an impact in nearly 80% of experiments, while the mature team will do so only in 25% of them!

What about win rates?

Show code 
estimate_win_rate = function(baseline, required_sample_size, density_func) {

  integrate(function(x) {

    p = power.prop.test(

      n=required_sample_size, 

      p1=baseline,

      p2=baseline + x

    )$power    

    density_func(x) * replace_na(p, 0)    

  }, 0, Inf)$value  

}



tibble(

  `Assumed effect size distribution` = sapply(eff_size_distrs, function(x) x$name),

  `P(detect a win | effect size distribution, N, baseline)` = sapply(

    eff_size_distrs,

    function(x) estimate_win_rate(baseline, required_sample_size, x$density)

  )

) |> kbl(caption='Long-run win rates', digits=2) |> kable_styling()
Long-run win rates
Assumed effect size distribution P(detect a win | effect size distribution, N, baseline)
Startup team 0.62
Mature team 0.23

We get 62% and 24%, respectively. The latter number isn’t far off an often-quoted figure for the win rate in online experiments of 10−20%10-20\%. That’s worth pondering: we know that the mature team produces winning ideas 75% of the time, but because most effects are minor, the actual observable win rate is just 24%.

Remember, these results were obtained using a sample size that “provides 80% power under an MDE of 1.5%”! But our “actual power” is closer to 30% (0.24/0.75). The MDE choice in the case of the mature team is clearly questionable.

Impact of sample size

How do win and detection rates change with the sample size? We can plot win and detection rates as a function of nn.

Show code 
sample_sizes = seq(1000, 200000, 2000)



win_and_detection_rates = pblapply(sample_sizes, function(n) {

  lapply(eff_size_distrs, function(d) {

    list(

      name = d$name, 

      detection_rate = estimate_detection_rate(baseline, n, d$density),

      win_rate = estimate_win_rate(baseline, n, d$density), 

      n=n

    )

  })

}, cl=4) |> bind_rows()



mde_sample_sizes = lapply(c(0.004, 0.005, 0.0075, 0.01, 0.015), function(mde) {

  list(

    n=power.prop.test(

      p1=baseline, p2=baseline + mde, power=0.8, sig.level=0.05

    )$n, 

    mde=mde

  )

}) |> bind_rows()



ggplot(win_and_detection_rates) + 

  geom_line(aes(x= n, y=detection_rate, color=name, linetype="Detection rate")) +

  geom_line(aes(x= n, y=win_rate, color=name, linetype='Win rate')) +

  geom_hline(aes(yintercept=0.75, color='Startup team', linetype='Share of positive-impact experiments')) +

  geom_hline(aes(yintercept=0.74, color='Mature team', linetype='Share of positive-impact experiments')) +

  theme_light() + 

  theme(

    panel.grid.major = element_blank(),

    legend.position = "right"

  ) + labs(

    x = 'Sample size (per group)',

    y = '',

    color = 'Effect size prior assumption',

    title = 'Avg detection and win rates as function of sample size',

    linetype = 'Metric'

  ) +

  scale_x_continuous(

    sec.axis = dup_axis(

      name = 'MDE that results in 80% power',

      breaks = mde_sample_sizes[['n']],

      labels = scales::percent(mde_sample_sizes[['mde']]),

      guide = guide_axis(n.dodge = 3)

    ),

    breaks = c(1000, seq(0, 200000, 5000)),

    labels = c("1,000", sapply(

      seq(0, 200000, 5000), 

      function(x) if_else(x %% 40000 == 0 && x > 0, format(x, big.mark = ","), ""))

    )

  ) +

  scale_y_continuous(breaks = seq(0, 1, 0.1), labels = scales::percent_format()) +

  scale_linetype_manual(

    values = c(

      "Detection rate"="solid", 

      "Win rate" = 'longdash',

      "Share of positive-impact experiments" = 'dashed'

    )

  )

There’s a lot going on in this chart, but here’s a few key takeaways:

  • The start-up team reaches ~90% detection rates (and, by extension, approaches its true win rate) as soon as the experiment sample size exceeds 40,000 observations per group. There’s little to be gained beyond that.

  • At the same time, if this were an early-stage start-up with just a few thousand customers to experiment with, it’s very possible it would see win/detection rates of just 50-60%. It’s a tough trade-off between the speed in decision-making that’s key in a start-up context and a higher share of false positives.

  • In the case of the mature team, there’s most to be gained until sample sizes reach ~80,000 observations per group. That’s an MDE of ~0.5%! Previously, we saw that using a 1.5% MDE (or 10,000 users per group) yields an observable win rate of 24%; increasing the sample size to 80,000 would double the observable win rate. Granted, most of these wins would be small, but as discussed above, that may still be a lot of money.

  • The mature team would have difficulty seeing win rates close to its true success rate (74%). It would take 5,000,000 observations per group to reach a win rate of 70%.

What about the Type I error rate in this context?

We saw how to go from statistical power / Type II error rate to long-term detection/win rates. Can we do the same with the Type I error rate?

Well, the thing is that the Type I error rate is only relevant if some experiments have zero impact. It’s not “close enough” to zero, but literally zero. In reality, (hopefully!) you don’t run such experiments. It works for individual experiments because the frequentist hypothesis testing framework asks “how likely we would have observed data like the one at hand if the true effect was zero” and not “how likely the true effect is zero.” But this question doesn’t translate to a series of experiments.

Additionally, the Type I error rate does not depend on the data you collect (sample size, variance). Even if you assume that some experiments have precisely zero impact, all you can do is choose Type I error rate with α\alpha; you can’t influence it beyond that choice.

As a result, there’s nothing we can optimize ¯\_(ツ)_/¯. You could ask deeper questions about the usefulness of an independent-of-data error rate that doesn’t apply in most situations, which would land you into the world of “Is null hypothesis testing the right thing to do?” but we’re not going there today.

What are the alternatives?

For one, we can estimate expected confidence interval widths and quantify the average uncertainty in our results. In a simple t-test setting: CIwidth=2MDE(β=0.5)CI_{\text{width}} = 2 MDE(\beta=0.5), i.e., the widths are twice the MDE that yields 50% power. Or, graphically:

Show code 
sample_sizes = seq(1000, 200000, 1000)



ci_widths = lapply(sample_sizes, function(n) {

  list(

    n=n,

    ci_width = (power.prop.test(p1=baseline, power=0.5, n=n)$p2 - baseline) * 2

  )

}) |> bind_rows()



ggplot(ci_widths) + 

  geom_linerange(aes(x=n, ymin=-ci_width/2, ymax=ci_width/2)) +

  theme_light() + 

  theme(

    panel.grid.major = element_blank(),

    legend.position = "right"

  ) + labs(

    x = 'Sample size (per group)',

    y = 'CI width',

    title = 'Confidence interval widths as a function of sample size'

  )

But in the spirit of focusing on “error rates,” let’s explore some other error types - meet Type S and Type M error rates.

Type M and S errors

Gelman and Carlin (2014) introduced two other metrics to use in statistical inference:

  • Type S (sign) error. It’s the probability that an observed difference, when declared statistically significant, has an opposite sign than the actual difference. In other words, it is the probability that you will make a decision that is opposite to what you are after. It directly relates to the practical implications of making a wrong decision. In a lot of situations, rolling out a zero-impact change is not the end of the world. Rolling out something that has an impact opposite to the one estimated in an experiment? That’s bad!
  • Type M (Magnitude) error measures how much the observed difference, when declared statistically significant, differs from the actual difference, expressed as a ratio. It partially addresses the issue of quantifying the risk of saying “effect sizes so small that they don’t practically matter” as statistically significant and is closely related to the winner’s curse phenomenon. I like Type M error a lot because it provides a rigorous quantification of the hand-wavy statement, “Yeah, the results are statistically significant, but I would not trust them due to the small sample size.” It provides a bridge between sample size and Type I error rate.

Calculating Type S and M error rates

The original paper proposed calculating Type S and M error rates via simulation. using a function. In 2019, Lu, Qui, and Deng published a paper that included closed-form formulas for these error rates. Their implementation is below (and can be easily extended to unequal sample size designs).

The rates are a function of four parameters: AA, the effect size; σ\sigma, the observed standard deviation of a metric; α\alpha, the chosen level of statistical significance; and nn, the sample size.

AA, conceptually, is a hypothesized effect size. That’s not the same as MDEMDE - arguably, in a well-designed experiment, MDEMDE should always be smaller than AA, as it’s the minimum effect size you’re interested in detecting.

Show code 
calculate_type_SM = function(A, sd, n, alpha = 0.05) {

  # A - effect size

  # sd - baseline standard deviation

  # n - size per group

  lambda = abs(A) / sqrt(sd**2/n + sd**2/n)

  z = qnorm(1 - alpha/2)

  neg = pnorm(-z -lambda)

  diff = pnorm(z - lambda)

  pos = pnorm(z + lambda)

  inv_diff = pnorm(lambda - z)

  list(

    power = neg + 1 - diff,

    typeS = neg / (neg + 1 - diff),

    typeM = (

      dnorm(lambda + z) + dnorm(lambda - z) + lambda * (pos + inv_diff - 1)

    ) / (lambda * (1 - pos + inv_diff))

  )

}

To get a sense of how these error rates look like, let’s use the setup from the examples above: a baseline of 17%, the hypothesized effect size of 1 percentage point, and a sample size of 10,000 users in each group:

Show code 
closed_form = calculate_type_SM(

  A=0.01, 

  sd=sqrt(baseline * (1 - baseline)), 

  n = required_sample_size

) |> as_tibble() |> mutate(method = 'Closed-form formula')



closed_form |>

  kbl(caption = 'Estimated Metrics', digits=6) |> kable_styling()
Estimated Metrics
power typeS typeM method
0.469165 0.00013 1.45038 Closed-form formula

The risk of a sign error in such a setting is very small - just ∼1/7700\sim1/7700. The type M ratio, on the other hand, is 1.451.45; a 45% overestimation can be meaningful in decision-making.

We can verify that these formulas are correct with a small simulation. In this simulation, we repeat the same experiment with the above parameters many times and estimate empirical error rate estimates. Note that estimating Type S error accurately via simulation isn’t easy given how rare it is, which is why even with 100,000 loops, we are still a bit off.

Show code 
pboptions(type="none")

sim_results = pbsapply(1:100000, function(i) {

  set.seed(i*24)

  A = rbinom(required_sample_size, 1, baseline)

  B = rbinom(required_sample_size, 1, baseline + 0.01)

  mean_diff = mean(B) - mean(A)

  p_value = t.test(B, A)$p.value  

  c(

    detected = p_value <= 0.05,

    sign_error = ifelse(p_value <= 0.05, sign(mean_diff) != sign(0.01), NA),

    ratio = ifelse(p_value <= 0.05, mean_diff / 0.01, NA)

  )  

}, cl=4) |> t() |> as_tibble()





empirical_rates = sim_results |> summarize(

  power = mean(detected, na.rm=T),

  typeS = mean(sign_error, na.rm=T),

  typeM = mean(ratio, na.rm=T),

  method = 'Empirical results'

)



bind_rows(closed_form, empirical_rates) |> 

  kbl(caption = 'Estimated Metrics') |> kable_styling()
Estimated Metrics
power typeS typeM method
0.4691649 0.0001298 1.450380 Closed-form formula
0.4581100 0.0000655 1.462017 Empirical results

Lu, Qiu, and Deng’s paper includes charts that perfectly illustrate the risk associated with underpowered experiments. On the x-axis, we plot statistical power (or, equivalently, sample size given a chosen A=MDEA=MDE), and on the y-axis - Type S error rate / Type M ratio:

Show code 
error_rates = lapply(seq(500, 40000, 500), function(n) {

  calculate_type_SM(0.01, sqrt(baseline * (1 - baseline)), n)

}) |> bind_rows() |> pivot_longer(-power)



ggplot(error_rates, aes(x=power)) + 

  geom_line(aes(y=value, color=name)) +

  scale_x_continuous(breaks=seq(0.05, 0.95, 0.2)) +

  facet_wrap(

    ~name, 

    scales='free_y', strip.position = 'top',

    labeller = as_labeller(c(typeM = 'Type M error ratio', typeS='Type S error rate'))

  ) +

  labs(x = 'Statistical power') +

  theme_light() + 

  theme(

    panel.grid.major.x = element_blank(),

    legend.position = "none",

    strip.placement = "outside"

  ) + ylab(NULL)

We can see that Type S error (which is very costly, decision-making-wise!) largely disappears once we hit the statistical power of 0.35−0.4. On the other hand, the type M ratio only goes below 1.5x once we achieve the statistical power of 0.50.5.

These charts should be a part of any Experimentation 101 course.

Obtaining long-term sign-error / exaggeration rates

Type M and S rates depend on a hypothesized effect size AA, just like the Type II error rate depends on MDEMDE, so we’ll need the same integration approach to get to long-term “average” rates. We’ll also tweak the definitions a bit and calculate:

  • Expected “sign error rate” - the % of experiments where we will declare the results statistically significant, but the observed treatment effect will have an opposite sign than the actual treatment effect. We’ll use all experiments as the denominator here, which is slightly different from the Type S error rate definition. If needed, we can divide this rate by the average detection rate to get back to the original definition.

  • Expected “average absolute exaggeration” - the expected absolute difference between the observed and actual treatment effects among experiments that we declare statistically significant. We’ll focus on absolute difference rather than a Type M-style ratio metric because integrating Type M ratio directly at small effect sizes is challenging as the ratio’s denominator approaches zero.

Show code 
estimate_sign_error_rate = function(baseline, s, density_func) {  

  baseline_sd = sqrt(baseline * (1 - baseline))

  integrate(function(x) {

    r = calculate_type_SM(n=s, A=x, sd=baseline_sd)

    (r$typeS * r$power) * density_func(x)    

  }, -Inf, Inf, stop.on.error = F)$value

}



estimate_abs_exaggeration = function(baseline, s, density_func) {  

  baseline_sd = sqrt(baseline * (1 - baseline))

  integrate(function(x) {

    r = calculate_type_SM(n=s, A=x, sd=baseline_sd)

    (r$typeM - 1) * density_func(x) * abs(x) * r$power    

  }, -Inf, Inf, stop.on.error = F)$value

}





sign_and_ratio_rates = pblapply(sample_sizes, function(n) {

  lapply(eff_size_distrs, function(d) {

    list(

      name = d$name, 

      sign_error_rate = estimate_sign_error_rate(baseline, n, d$density),

      exaggeration_ratio = estimate_abs_exaggeration(baseline, n, d$density), 

      n=n

    )

  })

  

}, cl=4) |> bind_rows() |> pivot_longer(

  cols=c('sign_error_rate', 'exaggeration_ratio'), names_to = "metric"

)

How do these rates look under the assumptions of effect size distributions we used earlier?

Show code 
ggplot(sign_and_ratio_rates) + 

  geom_line(aes(x= n, y=value, color=name)) +

  theme_light() + 

  theme(

    panel.grid.major = element_blank(),

    legend.position = "bottom"

  ) + labs(

    x = 'Sample size (per group)',

    y = '',

    color = 'Effect size prior assumption',

    title = 'Avg sign error rates and effect exaggeration size function of sample size',

    subtitle='Note: Discontinuities in the charts are an artifact of numeric integration...'

  ) +

  facet_wrap(

    ~metric, 

    scales='free_y', strip.position = 'top',

    labeller = as_labeller(c(

      exaggeration_ratio = 'Avg. absolute exaggeration', 

      sign_error_rate='Avg. sign error rate'

    ))

  ) +

  scale_x_continuous(

    breaks = c(1000, seq(0, 200000, 5000)),

    labels = c("1,000", sapply(

      seq(0, 200000, 5000), 

      function(x) if_else(x %% 40000 == 0 && x > 0, format(x, big.mark = ","), ""))

    )

  )

  • Broadly, once you enter 5000+ observations per group territory, sign error rates and average absolute exaggeration become limited. Even at “small” sample sizes (e.g., 1,000 observations), they are not worrisome, especially in the case of the start-up team.

  • The mature team has it worse on both metrics, given the same sample size; its experiments have a higher risk of sign error and a higher exaggeration ratio. This reflects the fact that, when the results are statistically significant, the effect sizes tend to be smaller (and thus at a higher risk of being influenced by a higher sampling error).

  • Again, a sample size of 80,000 observations per group seems to be where the mature team starts to see diminishing returns. The average absolute exaggeration decreases below 0.00050.0005 (5 basis points), and the sign error rate goes below 0.25%0.25\%.

Conclusion

So, what’s the TLDR?

  • A low observable win rate of (e.g., 25%25\%) doesn’t necessarily mean that the team produces winning ideas infrequently. It may simply be that the impact those ideas have tends to be smaller than you can reliably detect.

  • The only way to increase observable win/detection rates is by increasing sample sizes (or reducing the variance of the outcome metric). Choosing a higher MDE will not. If the effect sizes you tend to achieve are ~1ppt, you are artificially reducing your long-term win rate 2x if you run experiments with MDE of 1ppt instead of 0.5ppt.

  • Once you are in the 5,000+ sample size territory, the risk of sign errors is minimal. On the other hand, exaggeration ratios tend to remain sizable and should be considered when making decisions.

Does this mean a team operating in a mature setting should insist on 100,000+ observations in every experiment? It depends on how you answer questions like:

  • Do small effects truly matter? Perhaps having a 50% chance of detecting a 0.5ppt impact is OK? Maybe a tiny positive effect size is still a failure from a business perspective?

  • Do experiment runtimes impede the team’s ability to learn and execute? It’s one thing if you’re just impatient and need to accept that waiting for another four weeks will enable a more precise decision, bringing in $$$\$\$\$ for years to come. It’s another if you’re blocked by the experiment or run so few experiments that you’re not learning enough.

Most importantly, remember that fiddling with an MDE gives you nothing. Sample size is king. Use CUPED to reduce variance and sequential testing to enable peeking; otherwise, it’s all about the actual changes you test.

And, if you want to play around with effect size distrubution assumptions and see corresponding win/detection/sign error rates - I made a small Streamlit app just for that.

Appendix - Don’t trust the maths

Integration formulas look fancy, but do they work? I ran another simulation to double-check—it looks like it all adds up, with exception of exaggeration ratios that are a bit off. I suspect that differences arise from numerical approximations in integration.

Show code 
num_experiments = 100000



used_sample_size = floor(power.prop.test(

  p1=baseline, 

  p2=baseline + 0.015,

  power=0.2)$n)



appendix_sim = pblapply(eff_size_distrs, function(d) {

  drawn_effect_sizes = sapply(

    1:num_experiments, 

    function(i) d$generator()

  )  

  set.seed(42)

  exp_results = tibble(

    eff_size = drawn_effect_sizes,

    converted_A = rbinom(

      num_experiments, 

      used_sample_size, 

      baseline

    ),

    converted_B = rbinom(

      num_experiments, 

      used_sample_size, 

      baseline + drawn_effect_sizes

    )

  )|> mutate(

    conv_rate_A = converted_A / used_sample_size,

    conv_rate_B = converted_B / used_sample_size,

  ) |> mutate(

    mean_diff = conv_rate_B - conv_rate_A,

    st_error = sqrt(

      (conv_rate_A * (1 - conv_rate_A) / used_sample_size) +

      (conv_rate_B * (1 - conv_rate_B) / used_sample_size)

    )

  ) |>

    mutate(

      is_stat_sig = abs(mean_diff / st_error) >= qnorm(0.975)

    ) |> mutate(

      sign_error = if_else(

        is_stat_sig, 

        sign(mean_diff) != sign(drawn_effect_sizes), 

        NA

      ),

      exaggeration = if_else(is_stat_sig, abs(mean_diff - drawn_effect_sizes), NA)

    )

  

  exp_results |> summarize(

    distribution = d$name,

    detection_rate = mean(is_stat_sig),

    avg_sign_error = mean(replace_na(sign_error, F)),

    avg_exaggeration = mean(replace_na(exaggeration, 0)),

    avg_TypeS = mean(sign_error, na.rm=T),

    avg_conditional_exaggeration = mean(exaggeration, na.rm=T)

  )

}) |> bind_rows()



appendix_sim |> 

  kbl(caption='Simulation results using N = 20% power at MDE of 0.015', digits=4) |> 

  kable_styling()
Simulation results using N = 20% power at MDE of 0.015
distribution detection_rate avg_sign_error avg_exaggeration avg_TypeS avg_conditional_exaggeration
Startup team 0.5033 0.0028 0.0057 0.0056 0.0114
Mature team 0.0994 0.0125 0.0022 0.1262 0.0225
Show code 
lapply(eff_size_distrs, function(d) {

  list(

    distribution = d$name,

    detection_rate = estimate_detection_rate(

      baseline = baseline, 

      required_sample_size = used_sample_size, 

      density_func = d$density

    ),

    avg_sign_error = estimate_sign_error_rate(

      baseline = baseline, 

      s = used_sample_size, 

      density_func = d$density

    ),

    avg_exaggeration = estimate_abs_exaggeration(

      baseline = baseline, 

      s = used_sample_size, 

      density_func = d$density

    )

  )

}) |> bind_rows() |> 

  mutate(

    avg_TypeS = avg_sign_error / detection_rate,

    avg_conditional_exaggeration = avg_exaggeration / detection_rate

  ) |>

  kbl(caption='Integration results with equivalent parameters', digits=4) |> 

  kable_styling()
Integration results with equivalent parameters
distribution detection_rate avg_sign_error avg_exaggeration avg_TypeS avg_conditional_exaggeration
Startup team 0.4998 0.0026 0.0029 0.0052 0.0057
Mature team 0.0864 0.0118 0.0021 0.1364 0.0247
]]> 615 Estimating home court advantage in Lithuanian Basketball League with Gaussian Processes https://aurimas.eu/blog/2023/06/estimating-home-court-advantage-in-lithuanian-basketball-league-with-gaussian-processes/ Thu, 01 Jun 2023 04:04:05 +0000 https://aurimas.eu/blog/?p=602 notebook ]]> 602 Modeling tenure effects the Bayesian way https://aurimas.eu/blog/2023/04/modeling-tenure-effects-the-bayesian-way/ Thu, 27 Apr 2023 02:37:47 +0000 https://aurimas.eu/blog/?p=589
notebook
]]> 589 Getting faster to decisions in A/B tests – part 2: misinterpretations and practical challenges of classical hypothesis testing https://aurimas.eu/blog/2023/02/getting-faster-to-decisions-in-a-b-tests-part-2-misinterpretations-and-practical-challenges-of-classical-hypothesis-testing/ https://aurimas.eu/blog/2023/02/getting-faster-to-decisions-in-a-b-tests-part-2-misinterpretations-and-practical-challenges-of-classical-hypothesis-testing/#comments Thu, 16 Feb 2023 04:31:31 +0000 https://aurimas.eu/blog/?p=549

(this post is dedicated to Qixuan: thank you for all the patience on all the evenings I kept on bringing up various maddening questions, for the feedback on early drafts, and for all suggestions. This post wouldn’t be nearly as good if not for you!)

In the first post of the series, I covered some key approaches the industry uses to get to decisions faster than a classical null hypothesis testing approach allows. I also alluded that null hypothesis testing, while often highlighted for its basic statistical shortcomings (e.g., no peeking), has other pitfalls that are not immediately obvious.

In this post, I want to detail what I meant. In particular, I want to focus on practical challenges and easy-to-misinterpret concepts that are core to null hypothesis testing. I will not cover fallacies of p-values or issues in interpreting results – while they are confusing (here’s an academic paper that lists 25 (!) common misinterpretations), a lot of them come down to being pedantic and precise in communication. For the folks in the industry, there are more practical issues to consider.

A refresher on null hypothesis testing

Let’s begin with a quick recap of assessments of experiments using null hypothesis testing (NHT). Here is the process in a nutshell:

  1. Before the test begins, we must decide how long an experiment will run. NHT does not allow multiple evaluations of ongoing results, so we have to set a cut-off time. Why exactly can’t we peek? Some great blog posts explain the reasons already, but it comes down to the fact that peeking reduces specificity (increases the risk of false positives). Even a few peeks can increase error rates to 3-4x the intended rate. I discuss some statistical solutions to this issue in the first post of the series.
  2. The proper way to do it is by estimating the sample size required, which is a function of the certainty levels we want to achieve and the characteristics of the metric of interest. This is referred to as “power analysis”, and a bunch of calculators online help do the maths associating these factors, though it’s simple enough that you can easily code it up in a spreadsheet.
  3. Once we have the sample size, we know the time it will take to run the experiment, and the experiment can begin.

Four factors express the certainty levels and the characteristics of the metric of interest:

  • Desired sensitivity level. This is also known as power or 1 – Type II error rate (1 - \beta). It represents the true positive rate, i.e., the probability that the test will declare a difference between two variants when an actual difference exists between them. Usually, it is set to 1 - \beta = 0.8. The higher the sensitivity rate, the more data is required.
  • Desired specificity level, which represents the true negative rate. Typically, this is expressed as 1 – specificity, i.e., the error rate at which you declare variants in your experiment as different when they are the same. This is commonly referred to as the Type I error rate, a.k.a. statistical significance level \alpha. The most common values used are \alpha=0.05, 0.1. The higher the specificity level (the lower the error rate), the more data is required.
  • Minimum detectable effect (MDE) – the true effect size (difference in underlying populations) at which we want to achieve the specificity and sensitivity levels. Detecting a difference of 10% requires less data than detecting a difference of 1%. MDE typically depends on the expected impact of the change tested.
  • The variance of the metric. The more natural noise in the metric, the more data we need to differentiate the signal arising from the differences in the tested variants vs. the noise. In practice, variance is estimated using historical data of the metric of interest.

Now, let’s get into the challenges and misconceptions.

#1: How to choose specificity and sensitivity levels?

Statistics classes typically teach to use \alpha=0.05 as a default. Similarly, $\beta$ is usually set to $0.2$. Most software uses these default values, too, and I bet that most A/B tests in the industry do, too.

That is problematic. Should every business be 4x more concerned with false positives than false negatives by default? Is every experiment equal, and should these values be set in stone? Aren’t there some experiments where declaring an incorrect winner has limited business impact, but missing a difference may be more important? If not from a pure “did it move a metric” perspective, then from a team morale perspective?

Why do we use \alpha=0.05, \beta = 0.2? It’s not that they are special. For an overview of how they came to be, I recommend this paper, but it comes down to the personal preferences of influential statisticians that, over time, became the standard. This is not a robust basis for decision-making.

On top, as I will show later, specificity and sensitivity may not even mean what you think they mean. Specifically, they do not represent long-term error rates that you can expect to make across multiple experiments – which is a very intuitive way to think about them.

Are there alternatives that do not require thinking about sensitivity and specificity levels? Yes and no.

A surface answer is that Bayesian methods don’t rely on these exact concepts. However, there is no way to get away from quantifying uncertainty, and even Bayesian methods require some decisions about sensitivity and specificity levels. However, the way they go about it may be more intuitive, partly because it becomes part of interpreting results rather than being a discussion before the experiment even begins.

#2: How to choose MDE?

Challenge #1 is conceptual; most teams will likely “solve” it using default values. This challenge is more practical, something that product managers and data analysts actually think about. How do you choose the right MDE?

MDE represents the hypothesized true underlying difference in populations. So, conceptually, one should pick an MDE corresponding to the expected impact of the tested feature. But how do you reason about it? How do you decide if it should be 1% or 0.8% for a given feature? What about 0.5%? You could, in theory, use the track record of past experiments and knowledge of a given feature to guide MDE selection, but it’s really hard. After all, we run A/B tests precisely because we do not have a good gut feeling about the impact of the feature. If we knew it, there would be no need to test.

It does not help that MDE is, in practice, the only adjustable parameter in sample size determination. If touching error rates is out of the question, then one way to get results faster is by assuming a higher MDE. To make matters worse, the relationship between MDE and the sample size is non-linear. For example, a test with \alpha=0.05, \beta=0.8, where the expected conversion rate is 50%, an MDE of 2% will require ~10,000 observations in each variant. To detect 1% MDE, the sample size increases 4x. That may be a difference between running an experiment for a week vs. a month! There is a huge incentive to “pick” a higher MDE.

Choosing MDE is hard even if you can resist that incentive and stick to MDE as it should be thought about: the hypothesized difference in populations.

In practice, however, MDE is often misinterpreted as “impact we will be able to detect”, which further compounds the issue. Partly, I blame the very unintuitive MDE definition, which, I think, deserves a separate heading.

#3: What does MDE even mean?

A situation in the first post in this series illustrates a common misconception of what MDE means:

A product manager asks: “soo.. we always use a power level of 80%, and we ran these past five experiments with ~10,000 users each, and that implied a minimum detectable effect of 1.5% on our key metric.. but we detected an effect of 0.8% on two of them. How come? I thought the point of minimum detectable effect was that it’s a minimum change we can detect, but we seem to be able to detect lower changes, too?”

(this actually happened to me in real life)

I would not blame the product manager! My perspective of what MDE represents changed a lot while writing this post, too (and resulted in most redrafting..). Even when looking at first search results from the web, it “feels like” MDE has to do with the results you observe:

The minimum detectable effect (MDE) is the effect size which, if it truly exists, can be detected with a given probability with a statistical test of a certain significance level.

Analytics-toolkit.com

Minimum detectable effect (MDE) is a calculation that estimates the smallest improvement you are willing to be able to detect. It determines how “sensitive” an experiment is.

Optimizely.com

I think one should be excused to conclude that MDE is the smallest difference between two samples that produce p-values low enough to reject the null hypothesis. Or in other words, “MDE is the size of observed effect I can expect to detect with chosen error rates and given the sample size required”.

Except it’s not.

MDE has (almost) nothing to do with the sample differences you will observe. It’s an input parameter of your assumed true population differences to sample size/power calculations. The output of those calculations guarantees that you will reject a null hypothesis of no differences based on observed samples at selected error rates. Not a null hypothesis of differences larger than MDE; just a difference of zero (that is the answer to the product manager’s question).

This is not to say that MDE is useless. It forces you to make an assumption about the difference in underlying populations – exactly what you should be concerned with, instead of just the observed sample data. But given MDE may be the only “tunable parameter”, it’s easy for discussions about MDE to become more about “what impact we could detect” and not “what underlying impact we hypothesize.”

This line of thinking is very wrong. MDE should not be tunable to achieve statistical power levels – it’s like putting a carriage in front of a horse. Assuming a higher MDE for a low-impact feature doesn’t increase power. You cannot trade off MDE against sample size. The only choice available is sample size vs. power, not sample size vs. MDE.

Wouldn’t it be nice if there was an alternative that does not use this subject-to-abuse concept? Again, some Bayesian approaches (specifically, the expected loss metric) provide alternatives.

Bonus question

Is this whole MDE business clear to you? How about a bonus question?

Suppose you observe a statistically significant difference of 0.2 (e.g., p-value 0.04) in a typical experimental setup using an MDE of 0.5. How do we report results?

On the one hand, we observe a difference of 0.2, a p-value below our threshold, so we reject the null. We do not accept any alternative hypothesis – we just reject the one that they are equal. On the other hand, our sample size was set up in a way that should give this result 80% of the time for one particular alternative hypothesis (differences are at least 0.5).

So.. what is the effect you should report to the CEO? 0.2, the observed one? Or 0.5, which may be the true underlying effect? Could you say “we observed a difference of 0.2, there’s a 5% chance the result is by chance, and an 80% chance that the true difference is at least 0.5”?

Could you?

🤔

This brings me to the next point.

(the answer is – no; you can’t really report based on MDE. The 80% is the probability of observing the result given the alternative hypothesis; it’s not the probability of the alternative hypothesis given the result)

#4: NHT results are presented as binary

This challenge is of no fault of the statistical procedure itself but an important one to discuss. Most of the time, A/B test results are presented as binary: they are either (statistically) significant or not. If they are, the effect size is added (“feature X increased conversion by 0.5%”).

This sort of presentation hides the additional information we get from NHT – the statistical uncertainty. You may say – but what about confidence intervals? Sure, they may be used by some teams. But, likely, even when they are, they are only used when the effect is deemed statistically significant, to begin with.

In effect, thus, A/B test outcomes are binary. We either get enough evidence to reject the null hypothesis (typically, at a very high standard – allowing only a 5% error rate), or we don’t. This approach makes sense in, for example, a new drug research setting – the drug either works or does not. You want to be certain that one or the other is true.

But in a typical A/B test setting? Would it not be better to qualify a statement about outcomes with a certainty element? “we are 40% sure this variant increases conversions” or, even better, “we are 30% sure that this variant increases conversions by at least 1ppt; there also is a 15% chance it reduces conversions”.

One may say: sure, you can do that with NHT. Just show wider confidence intervals; 70%, 80% ones! Indeed, one could do that. My argument is not about statistics. It’s about the fact that, in practice, NHT usage withholds critical information from decision-makers. These people are paid to make the right decisions, and they should be able to make these choices with more information than the “significant or not” binary outcome of a test.

There is, however, a statistical element to it, too. NHT outcomes cannot be directly interpreted as probabilities. The confidence interval definition is weird, too. And then there’s the MDE business discussed above. This is where Bayesian methods again take the spot – interpreting results is more intuitive.

We’re left with just one more thing to cover, but it will take a bit more work and code to get through. Feel free to take a break!

#5: α and β are not long-term error rates

I glanced over one important aspect of sensitivity & specificity definitions earlier. They represent rates of avoiding errors. But rates over what, exactly?

Well, NHT is a procedure from frequentist statistics, and frequentist statistics are built on the concept of using samples to make inferences about the entire population. Correspondingly, these error rates reflect the probability of making an inference mistake purely due to the sampling error. The probability that we will get “unlucky” with a random draw from a population.

It’s easiest to illustrate with some simulations. Here is the setup:

  • Two variants: A (control) and B (treatment)
  • Decision metric: conversion rate (mean of a binary outcome variable)
  • Current conversion rate: 42%
  • Error rates: \alpha=0.1, \beta=0.2
  • We are comfortable with the minimum detectable effect (MDE) of 2ppt (i.e., we want to detect when the conversion rate of variant B drops below 40% or increases above 44%).

Plugging the above parameters into a power calculator yields a sample size of ~7,600 observations in each variant. 

#setup
alpha = 0.1
power = 0.8
baseline_conversion_rate = 0.42
minimum_detectable_effect = 0.02
power_calcs = power.prop.test(
  p1=baseline_conversion_rate, 
  p2=baseline_conversion_rate + minimum_detectable_effect,
  sig.level = alpha,
  power = power
)

n = ceiling(power_calcs$n)
print(power_calcs) # n = 7575

Scenario 1: no differences

Let’s start with a situation where variant B is equivalent to variant A. Here, I draw 7,600 observations for each variant from the Bernoulli distribution with p = 42%. Or, if you prefer, I randomly toss a coin with the same probability of landing on heads (42%). Then perform a t.test() and note if I found a difference at 1-\alpha=0.9 confidence level. I repeat it 10,000 times (as if I had an opportunity to redo the same experiment over and over) and summarize the results below.

no_simulations = 10000
variant_generator = function(n, p) rbinom(n=n, size=1, p=p)

results = sapply(1:no_simulations, function(i) {
  variant_A = variant_generator(n, baseline_conversion_rate)
  variant_B = variant_generator(n, baseline_conversion_rate)
  t_test = t.test(variant_A, variant_B)
  c(pvalue = t_test$p.value, statistic = t_test$statistic)
})
Observed outcomes under null

In ~10% of the cases, despite the hypothetical coin having the same underlying chance to land on heads in both populations, I observe that a null hypothesis test declares one of the variants better. That is exactly what \alpha represents: the expected false positive rate when we repeatedly draw samples from two identical underlying populations, and any differences in samples are purely due to chance.

In the real world, of course, we only draw the samples once for a given experiment, and NHT mechanics rely on the central limit theorem to figure out how likely the particular pair of samples we are looking at is different due to “bad luck” alone. The famous p-value represents that chance. It is a measure of surprise.

There are a couple of other things to keep in mind about how hypothesis tests work under null:

  • We are guaranteed a 10% false positive rate independent of the sample size. Despite the current sample size chosen to balance \alpha, \beta and MDE, we would also observe 10% false positive rates with smaller sample sizes.
  • p-values only measure the level of surprise after observing the data at hand, assuming that the null hypothesis (no difference) is true. That’s it.

The chart below shows one way to test if you understand what p-values mean. I plot the p-value distribution observed in the 10,000 simulations ran earlier. The distribution is uniform. If you can explain why, then congrats! You passed the p-value challenge.

P-value distribution under null

Scenario 2: there is an actual difference

What about \beta? Let’s simulate a situation where variants differ, and differ exactly by 2%, our MDE (let’s make variant B better). 

alternative_results = sapply(1:no_simulations, function(i) {
  variant_A = variant_generator(n, baseline_conversion_rate)
  variant_B = variant_generator(n, baseline_conversion_rate + minimum_detectable_effect)
  t_test = t.test(variant_A, variant_B)
  c(pvalue = t_test$p.value, statistic = t_test$statistic)
})

We detected that variant B is better – but only 80% of the time, exactly what the selected power rate promised. Under these perfect conditions, when the variants differ exactly by 2%, we have a 20% false negative error rate. If we were to reduce the sample size or change by how much the variants differ, we would also see a change in false negative rates.

Outcomes when variants are different

So what’s the issue?

These error rates represent the probability that we make a mistake in a given experiment due to “bad luck” – and, in the case of \beta, only assuming we correctly guessed the true underlying difference (and set MDE to it). But what about error rates across experiments?

It’s easy to intuitively extend this further and conclude that \alpha, \beta also represents the long-term error rates we are likely to make across multiple experiments. These sort of long-term error rates is something we should care about, right?

There is a catch, though. For \alpha, \beta to represent long-term error rates, every single experiment should have an MDE set perfectly to the true underlying difference. That’s unrealistic.

Instead, let’s explore more plausible scenarios. Let’s start with the best case – suppose we get the right MDE, on average. Here is how I will simulate it:

  • A single random draw in the simulation now represents one distinct experiment. We’re no longer in the hypothetical world where we repeatedly observe different draws of the same experiment; now, we only get to see results once.
  • Same \alpha, \beta;
  • A fixed MDE of 2% (because we think that, on average, that’s the impact of a feature, and we don’t have better information to set it more accurately for every experiment)
  • Because all power analysis inputs are the same, we have the same sample size in every simulation round.
  • In every experiment, Variant A will have a static success probability (42%). In reality, the control group success rate would change as we incorporate new features following A/B test results, but I will ignore it for simplicity.
  • Variant B will have a different success probability in every experiment. I will draw the feature’s impact from a normal distribution with a mean of 0.02 (corresponding to MDE) and a standard deviation of 0.02 (we hope that new feature effects, on average, are net positive and effects tend to be small). The success probability will thus be the baseline (42%) plus the feature impact.
real_life_results = sapply(1:no_simulations, function(i) {
  variant_A = variant_generator(n, baseline_conversion_rate)
  rv_mde = round(rnorm(1, mean=minimum_detectable_effect, sd=minimum_detectable_effect),2)
  variant_B = variant_generator(n, baseline_conversion_rate + rv_mde)
  t_test = t.test(variant_A, variant_B)
  c(pvalue = t_test$p.value, statistic = t_test$statistic, real_lift = rv_mde)
})

Here’s what the results look like. I rounded the random variable to two decimal places to have situations where the effect is “zero” (i.e., any effect between [-0.5%; 0.5%]).

We can see that:

  • When the real effect is “zero”, we incorrectly declare one of the variants as a better one ~10% of the time (yay, that’s the \alpha, great!)
  • When B is better, we get it right ~80% ( that’s the \beta!) of the time, with the remaining cases mostly being inconclusive decisions
  • When A is better, we get it right ~55% of the time 🤔.

Why such results for variant A? Well, that’s because B is, on average, better than A by 0.02 (MDE) in the cases it is better, but the reverse does not apply to A. If we only look at situations when A is better, the average effect is ~0.015. If we were to plug this value into a power calculator, together with the chosen sample size and \alpha, we would find that the power we get is ~60% – not too far off from what we observe. So that’s the first takeaway – the long-term error rates are not symmetric.

What we saw, however, is the best possible scenario – only because of the specific mean and standard deviation of the “feature impact” distribution that I chose. When these values are different, the error rates change. To illustrate, here is the outcome of experiments under different choices of \mu \text{ and } \sigma.

I hope this illustrates that \beta is not a long-term error rate. If, for example, you are early in the experimental journey, and half of your experiments hurt metrics, you may be able to detect it only ~65% of the time, despite running every test with 80% power. On the other hand, if your team is experienced and comfortable experimenting with large changes that can swing metrics either way (and you correctly guess the average impact), then your long-term detection rates are above 80%.

What about \alpha rate? It is ~90% in every simulation, right? Well… kind of. The issue with this error rate is that, unlike previously, the two samples are never drawn from an identical distribution. The feature impact is very unlikely to be exactly 0. And under such conditions, what do we consider “equal populations”? I got around that problem by rounding, and as long as rounding is relatively precise (2+ decimal points, based on my explorations), this error rate indeed hovers around 90%. Move it up to 1 decimal point, and this error rate becomes wonky, too.

Can we avoid rounding? Intuitively, it would make sense to invoke MDE here. Shouldn’t it represent the “effect we care about” / “effect below which we are ok to consider things equal”? I wish it did. But, as discussed in the MDE section, it is not what it means.

In the end, the point is that NHT, as a statistical procedure, does not make any promises about error rates across experiments with different populations. And if you care about long-term error rates, you’re out of luck. Are the approaches that provide guarantees for long-term error rates? I understand that the expected loss approach (see Convoy’s post) tries to address it, but I am somewhat skeptical that the results could be mathematically proven. So I’ll be keen to take it for a spin using similar tests to the above in subsequent posts.

Summary

Did I think I would write nearly 4,000 words about the most basic statistical test out there? Not really. But here I am. Apologies.

Let’s recap. One annoying statistical limitation of NHT is the inability to peek at the results. It’s easily fixable, though, and there are a lot of different procedures that provide alternatives.

Other issues are more conceptual. I argue that many practitioners may use default error rates that may not be the best choices for their work and that choosing MDE is difficult, not least because its meaning is quite convoluted.

Finally, my other callouts are mostly about the fact that frequentist statistics don’t give any promises about long-term error rates nor allow quantifying outcomes in a probabilistic fashion, which may be desirable in a business setting. To be clear, it’s not the fault of the t-test or frequentist statistics. Rather, it is the users of the statistical procedure making incorrect interpretations of the parameters of this procedure. But it is also true that these procedures do not deliver the type of information needed in the decision-making contexts discussed in this post.

Are there solutions to every one of these issues? I am not sure. Some are neatly handled in various Bayesian methods. However, those solutions aren’t better because “Bayesian statistics is inherently better” (I hope to show that the outcomes the two paradigms produce are not that different in a future post!), but because they shift the focus to different ways to quantify outcomes (or because they impose additional assumptions). That’s why I am much more excited about them than pure statistical solutions, such as “always valid p-values”.

But, as the next step, the plan is to cover sequential testing. Bayesian methods will be the des(s)ert.

All code used in this blog post is available on GitHub.

]]> https://aurimas.eu/blog/2023/02/getting-faster-to-decisions-in-a-b-tests-part-2-misinterpretations-and-practical-challenges-of-classical-hypothesis-testing/feed/ 1 549