<jupyter_start><jupyter_text>Testing the t testOne of the great advantages of using simulation is that you can test theassertions your teachers make.For example, in the [permutation and t-test page](permutation_and_t_test), we asserted that the t-test is notvalid when the underlying distribution of the numbers is not close to the[normal distribution](https://en.wikipedia.org/wiki/Normal_distribution).We can investigate this claim by simulating numbers from the null (ideal)world, and seeing what results we get from the t-test.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>The particular variant of the t-test that we were looking at in the page abovewas the *independent sample* t test for groups with similar variance.  Similarvariance means that the distribution of the values in the first group isroughly equal to the distribution in the second group.For example, soon we will be testing again for a mean difference between thenumbers of mosquitoes attracted to each of the 25 volunteers who drank beer,and the equivalent numbers for each of the 18 volunteers who drank water.See [the datapage](https://github.com/odsti/datasets/tree/master/mosquito_beer) formore details on the dataset, and [the data license page](../data/LICENSE).For an equal variance test, we assume that the spread of the beer values isroughly equal to the spread of the water values, as measured by the *standarddeviation*, or, equivalently, the *variance*.  Remember the variance is thesquared standard deviation.We can pull together the code in [permutation and t-testpage](permutation_and_t_test) to implement our own t-test.<jupyter_code># Import the Scipy statistics routines.
import scipy.stats as sps
def t_test(group1, group2):
    """ Independent sample t value and one-tail upper-tail p value.
    """
    g1_mean = np.mean(group1)
    g2_mean = np.mean(group2)
    omd = g1_mean - g2_mean  # The observed mean difference.
    errors = np.append(group1 - g1_mean, group2 - g2_mean)
    g1_n = len(group1)  # Number of observations in group1
    g2_n = len(group2)  # Number of observations in group2
    df = g1_n + g2_n - 2  # The "degrees of freedom".
    estimated_sd = np.sqrt(np.sum(errors ** 2) / df)
    t_stat = omd / (estimated_sd * np.sqrt(1 / g1_n + 1 / g2_n))
    upper_tail_p = 1 - sps.t.cdf(t_stat, df)
    return [t_stat, upper_tail_p]<jupyter_output><empty_output><jupyter_text>The only new thing in the implementation above is the second-to-last line,where we are using a routine in Scipy to calculate the probability value fromthe t statistic; the details of this are not important for our purpose.First we go back to the logic of this p value, which is very similar to thelogic for permutation test p values:* Notice that the function calculates `omd = np.mean(group1) -  np.mean(group2)`. Call `omd` the *observed mean difference*.* Assume that we are in the null (ideal) world where the numbers from `group1`  and the numbers from `group2` have been drawn at random from the *same*  distribution.* The p value is the probability, in this null world, of seeing a mean  difference that is equal to or greater than the observed mean difference  `omd`.You can also think of a p value as an *index of surprise*.  The p value tellsyou how often you would expect to see an observed mean different this large, orlarger, in the null (ideal) world.  If the p value is small, then the observedmean difference is surprising.  For example, if the p value is 0.05, it meansthat such difference only occurs 5% of the time by chance in the null world, or1 in 20 times.  You could say it was surprising at a 5% level.  Similarly a pvalue of 0.01 means the result would only occur 1% of the time in the nullworld, and it is surprising at a 1% level.Here we recreate the mosquito, beer, water data from the [permutation andt-test page](permutation_and_t_test):<jupyter_code>beer_activated = np.array([14, 33, 27, 11, 12, 27, 26,
                           25, 27, 27, 22, 36, 37,  3,
                           23,  7, 25, 17, 36, 31, 30,
                           22, 20, 29, 23])
water_activated = np.array([33, 23, 23, 13, 24,  8,  4,
                            21, 24, 21, 26, 27, 22, 21,
                            25, 20,  7, 3])<jupyter_output><empty_output><jupyter_text>We run our t-test over these data to get the same result you saw in the[permtuation and t-test page](permutation_and_t_test).<jupyter_code>t, p = t_test(beer_activated, water_activated)
print('t statistic:', t)
print('Upper-tail p value:', p)<jupyter_output><empty_output><jupyter_text>To check our function is doing the correct calculation, we show that the t andp values are the same as the ones we get from using the standard Scipy functionfor independent t-tests:<jupyter_code>result = sps.ttest_ind(beer_activated, water_activated)
print('Scipy t statistic:', result.statistic)
print('Scipy upper-tail p value:', result.pvalue / 2)<jupyter_output><empty_output><jupyter_text>Here is the observed difference in means:<jupyter_code># Observed mean difference
np.mean(beer_activated) - np.mean(water_activated)<jupyter_output><empty_output><jupyter_text>The t-test p value above asserts that a difference in means as large as theobserved difference, or larger, would only occur about 5% of the time in a null(ideal) world, where the beer and water values come from the same distribution.The observed result is surprising at around the 5% level.How would we check the assertion that the t-test is valid for normaldistributions?If it is valid, then consider the situation where we do in fact draw twosamples from *the same* normal distribution, and then ask the t test for a pvalue.  If the p value is 5%, it means that such a result should only occur bychance, in the null world, 5% of the time.So, we can repeat this procedure, drawing numbers that do in fact come from thenull world, and check that the t-test only tells us that the result issurprising at the 5% level --- about 5% of the time.<jupyter_code>n_iters = 10000
p_values = np.zeros(n_iters)  # Store the p values
for i in np.arange(n_iters):
    # Make 40 numbers from a normal distribution with mean 10, sd 2.
    # These are our numbers from the null world.
    randoms = np.random.normal(10, 2, size=40)
    # Split into two groups of size 20, and do a t-test.
    t, p = t_test(randoms[:20], randoms[20:])
    # Store the p value from the t-test.
    p_values[i] = p
# Show the first 5 p values.
p_values[:5]<jupyter_output><empty_output><jupyter_text>If the t-test calculation is correct, then we should only see a p value of 0.05or smaller about 5% of the time.<jupyter_code># Proportion of times the t-test said: surprising at 5% level.
np.count_nonzero(p_values <= 0.05) / n_iters<jupyter_output><empty_output><jupyter_text>Here the t-test is doing a good job --- it labels the result as surprising, atthe 5% level, about 5% of the time.Now we ask - does it matter if the group sizes are unequal?  To test this, wedo the same calculation, but split the numbers from the null world into onegroup of 3 and another of 37:<jupyter_code># t-test working on unequal group sizes.
p_values = np.zeros(n_iters)  # Store the p values
for i in np.arange(n_iters):
    # Make 40 numbers from a normal distribution with mean 10, sd 2.
    randoms = np.random.normal(10, 2, size=40)
    # Split into two groups of size 3 and 37, and do a t-test.
    t, p = t_test(randoms[:3], randoms[3:])
    # Store the p value from the t-test.
    p_values[i] = p
# Show the first 5 p values.
p_values[:5]<jupyter_output><empty_output><jupyter_text>How good a job is it doing now, with unequal group sizes?<jupyter_code># Proportion of times the t-test said: surprising at 5% level.
# This time wih unequal group sizes.
np.count_nonzero(p_values <= 0.05) / n_iters<jupyter_output><empty_output><jupyter_text>The proportion is still around 5%, close to what it should be.What happens if we use a distribution other than the normal distribution?Here we use some random numbers from a [Chi-squareddistribution](https://en.wikipedia.org/wiki/Chi-squared_distribution).  Thedistribution looks like this, with a $k$ value of 2 (see the Wikipedia page):<jupyter_code>some_chi2_numbers = np.random.chisquare(2, size=1000)
plt.hist(some_chi2_numbers)
plt.title('1000 random samples from chi-squared distribution, k=2');<jupyter_output><empty_output><jupyter_text>We use this highly not-normal distribution to provide numbers to our t-test:<jupyter_code># t-test working on unequal group sizes and not-normal distribution.
p_values = np.zeros(n_iters)  # Store the p values
for i in np.arange(n_iters):
    # Make 40 numbers from a chi-squared distribution with k=2
    randoms = np.random.chisquare(2, size=40)
    # Split into two groups of size 3 and 37, and do a t-test.
    t, p = t_test(randoms[:3], randoms[3:])
    # Store the p value from the t-test.
    p_values[i] = p
# Show the first 5 p values.
p_values[:5]<jupyter_output><empty_output><jupyter_text>In this situation the t-test starts to be less accurate - labeling too manyrandom differences as being surprising at the 5% level:<jupyter_code># Proportion of times the t-test said: surprising at 5% level.
# This time wih unequal group sizes.
np.count_nonzero(p_values <= 0.05) / n_iters<jupyter_output><empty_output><jupyter_text>Does a permutation test do a better job in this situation?We can test!Here is a function that does a permutation test:<jupyter_code>def permutation(group1, group2, niters=10000):
    omd = np.mean(group1) - np.mean(group2)
    g1_n = len(group1)
    fake_mds = np.zeros(niters)
    pooled = np.append(group1, group2)
    for i in np.arange(niters):
        shuffled = np.random.permutation(pooled)
        fake_mds[i] = np.mean(shuffled[:g1_n]) - np.mean(shuffled[g1_n:])
    return np.count_nonzero(fake_mds >= omd) / niters<jupyter_output><empty_output><jupyter_text>Test this on the mosquito data:<jupyter_code>permutation(beer_activated, water_activated)<jupyter_output><empty_output><jupyter_text>This is very similar to the t-statistic p value --- *for these data* that havefairly equal group size, and a distribution not far from normal:<jupyter_code>t_test(beer_activated, water_activated)<jupyter_output><empty_output><jupyter_text>Now let's check how the permutation test does when there are unequal groupsizes and a not-normal distribution.The code below will take a few tens of seconds to run, because you are runningmany loops in the `permutation` function, each time you go through the mainloop.<jupyter_code># Permutation working on unequal group sizes and not-normal distribution.
# This is slow - do fewer iterations.
n_iters = 1000
p_values = np.zeros(n_iters)  # Store the p values
for i in np.arange(n_iters):
    # Make 40 numbers from a chi-squared distribution with k=2
    randoms = np.random.chisquare(2, size=40)
    # Split into two groups of size 3 and 37, and do a t-test.
    # Use fewer iterations than usual to save computation time.
    p = permutation(randoms[:3], randoms[3:], niters=1000)
    # Store the p value from the permutation test.
    p_values[i] = p
# Show the first 5 p values.
p_values[:5]<jupyter_output><empty_output><jupyter_text>How does the permutation test do?<jupyter_code># Proportion of times the permutation test said: surprising at 5% level.
# With unequal group sizes, not-normal distribution.
np.count_nonzero(p_values <= 0.05) / n_iters<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08075019536337588, 'fraction_numerical': 0.022922636103151862, 'mean_word_length': 4.3348772579898105, 'pattern_counts': {'":': 0, '<': 77, '<?xml version=': 0, '>': 74, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '407006', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3799, 'n_tokens_neox': 3514, 'n_words': 1467, 'path': 'permutation/testing_t.ipynb', 'repo_name': 'uob-ds/cfd2021', 'size': '11517'}