I want to estimate what sample size I need in an experiment where I compare the mean of an almost exponentially distributed variable between two groups (in an A/B test).

I know I'm always going to have more than 500 samples per variation so I assume that the difference in mean is normally distributed and calculate $z$ as

$$z = \frac{\bar {X}_1 - \bar{X}_2}{s_p \cdot \sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}$$

with $$ s_p = \sqrt{\frac{(n_1-1)s_{X_1}^2+(n_2-1)s_{X_2}^2}{n_1+n_2-2}}$$

To calculate the power I followed the Wikipedia example but replaced their standard error, $\hat{\sigma}_D/\sqrt{n}$, with $$SE=s_p\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}$$ which means that I end up with $$1-\beta = power\approx 1- \Phi(z_{\alpha/2}-\frac{\theta}{SE})$$ where $\Phi$ is the cumulative distribution function of the normal distribution and $\theta$ is the minimum difference in mean I want to be able to detect with rate $1-\beta$.

For simplicity I assume that $n_1=n_2$ and end up with the expression $n=\frac{2(s_{x_1}^2+s_{x_2}^2)(Z_{\alpha/2}-Z_{1-\beta})^2}{\theta^2}$

To validate this I did the script below where I generate three exponential distributions, one of them, A, with a scale parameter that is $mde$ (minimum detectable effect) higher compared to the other two distributions (B1 and B2). Then in each loop I do two "experiments" where I in one case compare two distributions with the same parameters and one where I compare two distributions with a small difference. Then when counting the number of "significant differences" in the experiments I expect to get false positive rate $\alpha$ when comparing B1 and B2 and true positive rate $1-\beta$ when comparing A with B1. When doing this on sampled distributions where I know all the parameters it seems to work. I get the power and significance I'm using as input.

Now the question is how I would apply this to a real scenario where I only can estimate $s_{x_1}$ and $s_{x_2}$? I tried to do this with some real data by using the standard deviation of all the data for both $s_{x_1}$ and $s_{x_2}$ and ended up really wrong in my sample size calculations (I got a way too high power). Any suggestions on how to approach this in a better way?

import numpy as np
import scipy.stats as st
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.optimize import fsolve, root
%matplotlib inline

def two_sample_z_test(sample1, sample2):
    s_p = np.sqrt(((sample1.size - 1) * sample1.std()**2 + (sample2.size - 1) * sample2.std()**2) /
    z_value = (sample1.mean()-sample2.mean())/(s_p * np.sqrt(1/sample1.size + 1/sample2.size))
    p_value = st.norm.sf(np.abs(z_value), sample1.size + sample2.size - 2) * 2  
    return p_value, z_value

def confidence_interval_normal(sample1, sample2,  alpha = 0.05):
    n1 = sample1.size
    n2 = sample2.size
    std1 = sample1.std()
    std2 = sample2.std()
    df = n1 + n2 - 2
    z_crit = st.norm.isf(q=alpha/2)
    s_p = np.sqrt(((n1 - 1) * std1**2 + (n2 - 1) * std2**2)/df)
    mean_diff = sample1.mean() - sample2.mean()

    ci_width = z_crit*s_p*np.sqrt(1/n1+1/n2)
    ci_lower = mean_diff - ci_width
    ci_upper = mean_diff + ci_width
    return ci_lower, ci_upper

def sample_size_z_test(mde, std1, std2, alpha=0.05, power=0.8):
    z_crit = st.norm.isf(q=alpha/2)
    z_power = st.norm.isf(q=power)

    n = 2*(std1**2+std2**2)*(z_crit-z_power)**2/mde**2

    return n

#Parameters for the distributions
b_mean = 10
mde = 1
a_mean = b_mean + mde
alpha = 0.05
power = 0.8

#Calculate power
sample_size=int(sample_size_z_test(mde=mde, std1=a_mean, std2=b_mean, alpha=alpha, power=power)/2)

n_loops = 5000 #number of "experiments" to run
ci_diff = np.zeros([n_loops,2])
ci_no_diff = np.zeros([n_loops,2])

diff_A_B1 = np.zeros(n_loops)
diff_B1_B2 = np.zeros(n_loops)

za = np.zeros(n_loops)
zb = np.zeros(n_loops)

for i in range(n_loops):
    experiment_A = np.random.exponential(size=sample_size, scale=a_mean)
    experiment_B1 = np.random.exponential(size=sample_size, scale=b_mean)
    experiment_B2 = np.random.exponential(size=sample_size, scale=b_mean)

    pa, za[i] = two_sample_z_test(experiment_A, experiment_B1)
    pb, zb[i] = two_sample_z_test(experiment_B1, experiment_B2)

    ci_diff[i,:] = confidence_interval_normal(sample1=experiment_A, sample2=experiment_B1, alpha=alpha)
    ci_no_diff[i,:] = confidence_interval_normal(sample1=experiment_B1, sample2=experiment_B2, alpha=alpha)

    diff_A_B1[i] = experiment_A.mean() - experiment_B1.mean()
    diff_B1_B2[i] = experiment_B1.mean() - experiment_B2.mean()

#Count how many times we are outside the 95% confidence interval (i.e. a significnace change=positive result)
significant_diff = (ci_diff[:,0]>0) | (ci_diff[:,1]<0)
significant_no_diff = (ci_no_diff[:,0]>0) | (ci_no_diff[:,1]<0)

# In how many "experiments" did we get a confidence interval that did not contain 0?
# This is the false positive rate for the case when there is no difference and true positive rate when 
# there is a difference
positive_rate_improvement = sum(significant_diff)/n_loops
positive_rate_no_diff = sum(significant_no_diff)/n_loops

print('True positive rate: ', positive_rate_improvement, '(chosen power: ', power, 
      ')\nFalse positive rate comparing B1 to B2: ', positive_rate_no_diff, 
      ' (chosen significance level: ', alpha,')')

ax1 = plt.subplot(1,1,1)
sns.distplot(diff_A_B1, ax=ax1, label = 'Difference in mean between A and B1')
sns.distplot(diff_B1_B2, ax=ax1, label = 'Difference in mean between B1 and B2')

x = np.linspace(-6,6,1000)
z_theory = st.norm.pdf(x=x)

ax2 = plt.subplot(1,1,1)

sns.distplot(zb, ax=ax2, label = 'actual z values')
ax2.plot(x, z_theory, label='theoretical z-curve')

EDIT It is an exaggeration to say that the data is exponentially distributed. It looks something like this: Distribution of data

  • $\begingroup$ Since the mean in a one-parameter exponential is a scale parameter, the more typical way to compare means of exponential distributions would be as a a ratio (does the ratio of scale parameters differ from 1?). While in sufficiently large samples you would indeed have $\bar{X}$ being approximately normal, the denominator in the t-test will not be independent of the numerator and variance estimate in the denominator will not be the variance times a chi-squared (so the statistic will not have a t-distribution); you'd require larger samples than needed for the numerator to be reasonably normal $\endgroup$ – Glen_b -Reinstate Monica Jul 24 '17 at 12:37
  • $\begingroup$ ... but even if you had a sample large enoug that you didn't have to worry about any of that, your power will be lower than some other tests (that is, for a given level of power your chosen test will require a larger sample than a more wisely chosen test). $\endgroup$ – Glen_b -Reinstate Monica Jul 24 '17 at 12:39
  • $\begingroup$ How do you know this variable is exponentially distributed? $\endgroup$ – Glen_b -Reinstate Monica Jul 24 '17 at 12:39
  • $\begingroup$ My bad, it is not really exponentially distributed (see new fig above). I don't know the true distribution but I thought it was OK to assume that the mean of the distribution will be normally distributed even though the data definitely is not? What would be a more wisely chosen test? Mann-Whitney U test? $\endgroup$ – Alex Jul 24 '17 at 13:25
  • $\begingroup$ In sufficiently large samples you will have that the difference in means would be close to normally distributed. Let's assume you have $n$ large enough for this to be the case. That's the numerator of the t-statistic (your "z"), but the t-statistic is not just a numerator; it's the behaviour of the ratio that matters. In really large samples you can ignore the distribution of the denominator and just treat the entire expression as normal (a Z test). However, this normality simply means that the test will have the right significance level. It will not necessarily mean you have the best power. $\endgroup$ – Glen_b -Reinstate Monica Jul 24 '17 at 17:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.