<jupyter_start><jupyter_text>Week 2 - Lab Session 3In this lab session we are going to look at how to calculate confidence intervals and run a t-test using Python. As per last week, you should start by running the existing cells which will load in some libraries and a required dataset. Note that these cells may throw a 'FutureWarning' which can be safely ignored.<jupyter_code>import statsmodels.api as sm
import pandas as pd
import scipy as sp
import numpy as np
df = sm.datasets.cpunish.load_pandas().data
df.describe()
df.head()<jupyter_output><empty_output><jupyter_text>Let's first get a mean and confidence interval for the income variable in the dataset. We can do that with the following set of commands:<jupyter_code>mean = df['INCOME'].mean()
se = df['INCOME'].std()/np.sqrt(df['INCOME'].shape[0])
sp.stats.norm.interval(0.95, loc=mean, scale=se)<jupyter_output><empty_output><jupyter_text>Check that you understand how this is calculated! Write some code which approximately duplicates the numbers above. Recall the formula for an approximate 95% confidence interval:\begin{equation*}95\%\ CI\ =\ mean\ ±\ (1.96\ *\ standard\ error)\end{equation*}<jupyter_code>(mean - (1.96 * se), mean + (1.96 * se))<jupyter_output><empty_output><jupyter_text>The amount of observations in this dataset is quite small, which means it would be important to use a t-distribution to calculate our confidence interval. We can do that with `sm.stats.DescrStatsW(df['INCOME']).tconfint_mean()`. How does using the t-distribution change the confidence interval?<jupyter_code>sm.stats.DescrStatsW(df['INCOME']).tconfint_mean()<jupyter_output><empty_output><jupyter_text>Now let's have a look at how to do a t-test. We're going to look at whether income is systematically higher or lower in southern states in the dataset. First, write a little bit of code that works out the difference in means between the two states. To get you started, you could use `df[df['SOUTH']==1]['INCOME'].mean()` to get the mean of all states which have south = 1.<jupyter_code>df[df['SOUTH']==1]['INCOME'].mean() - df[df['SOUTH']==0]['INCOME'].mean()<jupyter_output><empty_output><jupyter_text>Now let's run the t-test. Statsmodels has a library to do this for us: `sm.stats.ttest_ind(df[df['SOUTH']==1]['INCOME'], df[df['SOUTH']==0]['INCOME'], usevar="unequal")`. You can see the documentation here: https://www.statsmodels.org/dev/generated/statsmodels.stats.weightstats.ttest_ind.htmlThis returns three values: the test statistics, the p value, and the degrees of freedom<jupyter_code>sm.stats.ttest_ind(df[df['SOUTH']==1]['INCOME'], df[df['SOUTH']==0]['INCOME'], usevar="unequal")<jupyter_output><empty_output><jupyter_text>Challenge: write out a function that performs a t test and returns the test statistic. This will involve taking into two arrays of numbers, calculating the difference in means, and expressing that difference in means in terms of a number of standard errors. Check your result agrees with the function above. If you want to do even more, you could also calculate the degrees of freedom (see formula here: https://www.statisticshowto.datasciencecentral.com/satterthwaite-formula/) and hence determine the p value. Note that you can use `sp.stats.t.cdf(ts, df=deg_freedom)*2` to convert a test statistics (`ts`) and degrees of freedom (`deg_freedom`) into a p value.<jupyter_code>def calc_dof(x, y): 
    dof_num = (x.var()/x.shape[0] + y.var()/y.shape[0])**2 
    dof_denom = ((x.var()/x.shape[0])**2 / (x.shape[0]-1) + (y.var()/y.shape[0])**2 / (y.shape[0]-1))
    return(dof_num/dof_denom)
def t_test(df=None, var=None, group=None):
    #check there are two groups
    groups = df[group].unique()
    if len(groups)!=2:print('error - group variable must have exactly two values');return
    
    #calculate group means, variance and n
    means = [df[df[group]==g][var].mean() for g in groups]
    variances = [df[df[group]==g][var].var() for g in groups] 
    ns = [df[df[group]==g][var].shape[0] for g in groups]
    
    se_2 = (variances[0] / ns[0]) + (variances[1] / ns[1])
    se = np.sqrt(se_2)
    
    deg_freedom = calc_dof(*(df[df[group]==g][var] for g in groups))
    t = (means[0] - means[1]) / se
    p = sp.stats.t.cdf(t, df=deg_freedom)*2
    print('Comparing groups:', groups, 'for variable', group)
    print('Mean for group', groups[0], ':', means[0])
    print('Mean for group', groups[1], ':', means[1])
    print('Difference in means:', means[0] - means[1])
    print('Test statistic:', t)
    print('Degrees of freedom:', deg_freedom)
    print('P value:', p)
    
t_test(df=df, var='INCOME', group='SOUTH')<jupyter_output>Comparing groups: [1. 0.] for variable SOUTH
Mean for group 1.0 : 33325.42857142857
Mean for group 0.0 : 36958.8
Difference in means: -3633.3714285714304
Test statistic: -1.587305006418155
Degrees of freedom: 13.885455969119223
P value: 0.1349473046162541
{'fraction_non_alphanumeric': 0.13536137198856676, 'fraction_numerical': 0.03082890975908534, 'mean_word_length': 5.048148148148148, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '206977', 'license': 'MIT', 'n_tokens_mistral': 1726, 'n_tokens_neox': 1620, 'n_words': 525, 'path': 'Week_2_Lab_AA.ipynb', 'repo_name': 'abdulra4/Analytical-Statistics', 'size': '4898'}