# A/B Test from two distributions with unknown means

I will run a experiment to check if a given change, represent some significance in user interactions. I have two sets, sampled from two groups, one called control which counts the number of interactions from 30 sampled users in control group, and I have treatment which counts the number of interactions from 30 sampled users with a test feature. They are given in the Python code below

import numpy as np

control   = np.array([1, 0, 1, 3, 2, 1, 0, 1, 3, 2, 1, 0, 1, 3, 2, 1, 0, 1, 3, 2, 1, 0, 1, 3, 2, 1, 0, 1, 3, 2])
treatment = np.array([0, 1, 3, 2, 1, 2, 1, 3, 2, 1, 0, 2, 3, 2, 1, 0, 2, 3, 2, 1, 0, 2, 3, 2, 1, 0, 2, 3, 2, 4])


For this test i will use the following parameters

Practical Significant Boundary: $$d_{min} = 0.05$$

Significance level $$\alpha = 0.05$$

d_min = 0.05
alpha = 0.05


I want to check if in average there is no usage difference between the two groups. Then I state the null hypothesis as

$$H_0: \mu_t == \mu_c \implies d = 0$$

where $$\mu$$ is the average of control group (c) or treatment group (t) and $$d$$ the difference $$d = \mu_t - \mu_c$$.

mu_c = control.mean()
mu_t = treatment.mean()
d = mu_t - mu_c


By the characteristics of the data, I will use the t-test for evaluation. I will assume the variance is the same between both groups, that will give me the following computation to test statistics.

n_c = len(control) # Number of samples in control
n_t = len(treatment) # Number of samples in treatment

df = n_c - 1 + n_t - 1 # Degree of freedom

# squared sum of errors from control group
SS_c = np.sum((control - mu_c)**2)

# squared sum of errors from treatment group
SS_t = np.sum((treatment - mu_t)**2)

# Pooled standard error
S_pool = np.sqrt((SS_c + SS_t)/df)

# Test statistics for t-distribution
TS = d / (S_pool * np.sqrt(1/n_c + 1/n_t))


Now I will verify if the test confirms or rejects the null hypothesis. For that I compute the t_alpha with 58 degrees of freedom,

import scipy.stats as stats

t_alpha = stats.t.ppf(1-alpha/2, df)

TS > t_alpha or TS < -t_alpha
# Returns true


Since TS is out of the boundaries of t_alpha we can reject the null hypothesis.

Now I will verify if the result is practically significant, that is if my confidence interval including the standard error intersects somehow the significance region. For that I need to estimate the margin error ($$m$$), then compute the confidence interval from $$d \pm m$$

m = t_alpha * S_pool
ci = np.array([d - m, d + m])


Just for reminder d is the difference between the averages from two groups (control and treatment).

To visualize we can build the given plot.

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10, 3))
ax.vlines(0, -1, 1)
ax.vlines([-d_min, d_min], -1, 1, linestyles='dashed', label="Significant boundary");
ax.plot(d, 0,'b.', label="Center Confidence Interval");
ax.hlines(0, *ci, 'b', label="Confidence Interval");
ax.set_title("Significance");
plt.legend();


Since the confidence interval overlaps with the significant boundary (totally), even cross the $$0$$, I cannot assume that this test has practical significance, hence I would reject it.

The problem is: I was not expecting an interval SO large. Probably I just have used a wrong formula for margin error, but it seems correct. Any help to go forward in this is welcome.

• Why do you think it is wrong? It seems reasonable. You have very few observations. Try just using an online calculator to verify. Jan 14 at 4:55

It's always a good idea to use existing implementations to check, especially for something as common as a t-test.

import statsmodels.stats.weightstats as sm
test = sm.CompareMeans(sm.DescrStatsW(treatment), sm.DescrStatsW(control))
test.ttest_ind()
>>> TtestResult(statistic=1.0934, pvalue=0.2788, df=58.0)

test.tconfint_diff()
>>> (-0.2492, 0.8492)


Note that both contradict the results you have obtained. A statistical language like R will confirm this and give you everything you want to know in a single line of code:

t.test(treatment, control, var.equal = TRUE)
>   Two Sample t-test
> t = 1.0934, df = 58, p-value = 0.2788
> 95 percent confidence interval:
>  -0.2492  0.8492


Clearly you have at least one mistake in your formulae. I'm not sure why you end up with TS outside of the critical values because for me this is not the case when re-running your code:

TS = d / (S_pool * np.sqrt(1/n_c + 1/n_t))  # 1.0934
t_alpha = stats.t.ppf(1-alpha/2, df)        # 2.0017
TS > t_alpha or TS < -t_alpha
>>> False


The problem in your confidence interval calculation is that what you've called "pooled standard error" is actually the pooled standard deviation. The standard error is this quantity scaled by the square root of $$n$$ and also the denominator of your t statistic: S_pool * np.sqrt(1/n_c + 1/n_t). Using that you'll end up with the correct CI:

SE_pool = S_pool * np.sqrt(1/n_c + 1/n_t)
m = t_alpha * SE_pool
ci = np.array([d - m, d + m])
>>> array([-0.2492,  0.8492])


Finally, it's not clear what you mean by 'practical significance boundary'. The most liberal interpretation (and how you started your implementation) is that you want to test $$H_0: \mu_t-\mu_c = 0$$ and if you observe d > d_min you conclude that the effect might be practically relevant - if it were statistically significant. A more stringent interpretation, where you seem to want to go with your confidence interval, is that you actually want to test $$H_0: \mu_t-\mu_c\le d_{min}$$ which is the same as a $$1-\alpha$$ interval excluding d_min. I don't know which of these two you actually want, just highlighting that they are not the same. Neither of the null hypotheses can be rejected at the specified $$\alpha$$ in your case.

• Indeed. I made a small mistake while transcribing the data from my local file to the SO. And here, at StackOverflow, with the "wrong" data, the test statistics confirms the null hypothesis. Good point is that using both Python Libs and R the test statistics values (the TS variable), are the same. it means that equation is correct. And indeed, my misnaming took me to a conceptual error, thank you for catching this up. With that fix (as you shown) the results are now equal than R and Python Stats. BTW, is there a name for this $np.sqrt(1/n_c + 1/n_t)$? It looks like a kind of harmonic mean.
– Lin
Jan 15 at 5:09
• If you keep the variances (which you already pooled in S_pool) in the denominators they are a weighted pooled variance. It's usually this variance that you divide by $n$ (or $n-1$) and the standard error by the degrees of freedom and not the other way around - the end result is the same regardless of order though. Jan 15 at 7:36