I am comparing stats.ttest_ind() vs “manual” computation of the same test, and get different results.
import numpy as npimport pandas as pdimport scipy.stats as statsimport mathstats.ttest_ind() method:
#generate datanp.random.seed(123)df = pd.DataFrame({ 'age':np.random.normal(40,5,200).round(), 'sex':np.random.choice( ['male', 'female'], 200, p=[0.4, 0.6]), })#define groupsmen = df.age[df.sex == 'male']women = df.age[df.sex == 'female']#run t-testtest_stat, test_p = stats.ttest_ind(men, women)print(test_stat, test_p)Out:
-0.9265613940505325 0.355282312357339Manual method:
#meanmen_mean, women_mean = men.mean(), women.mean()#standard deviationmen_sd, women_sd = men.std(ddof=1), women.std(ddof=1)#standard errormen_n, women_n = len(men), len(women)men_se, women_se = men_sd/math.sqrt(men_n), women_sd/math.sqrt(women_n)#standard error on the difference between men and womense_diff = math.sqrt(men_se**2.0 + women_se**2.0)#t-statt_stat = (men_mean - women_mean) / se_diff#degrees of freedomdf = men_n + women_n - 2#critical valuealpha = 0.05cv = stats.t.ppf(1.0 - alpha, df)# p-valuep = (1 - stats.t.cdf(abs(t_stat), df)) * 2print(t_stat, cv, p)Out:
-0.9244538916746341 0.3563753194455255We can see there’s a small difference. Why? Maybe because of how stats.ttest_ind() computes degrees of freedom? Any insight much appreciated.
Advertisement
Answer
The following works. It is your code from above, with only two rows changed.
import numpy as npimport pandas as pdimport scipy.stats as statsimport math#generate datanp.random.seed(123)df = pd.DataFrame({ 'age':np.random.normal(40,5,200).round(), 'sex':np.random.choice( ['male', 'female'], 200, p=[0.4, 0.6]), })#define groupsmen = df.age[df.sex == 'male']women = df.age[df.sex == 'female']#run t-test############################### CHANGED THE ROW BELOW HEREtest_stat, test_p = stats.ttest_ind(men, women,equal_var=False) print(test_stat, test_p)#meanmen_mean, women_mean = men.mean(), women.mean()#standard deviationmen_sd, women_sd = men.std(ddof=1), women.std(ddof=1)#standard errormen_n, women_n = len(men), len(women)men_se, women_se = men_sd/math.sqrt(men_n), women_sd/math.sqrt(women_n)#standard error on the difference between men and womense_diff = math.sqrt(men_se**2.0 + women_se**2.0)#t-statt_stat = (men_mean - women_mean) / se_diff#degrees of freedom############################### CHANGED THE ROW BELOW HEREdf = (men_sd**2/men_n + women_sd**2/women_n)**2 / ( men_sd**4/men_n**2/(men_n-1) + women_sd**4/women_n**2/(women_n-1) )#critical valuealpha = 0.05cv = stats.t.ppf(1.0 - alpha, df)# p-valuep = (1 - stats.t.cdf(abs(t_stat), df)) * 2print(t_stat, cv, p)and it outputs
-0.9244538916746341 0.356441636045986-0.9244538916746341 1.6530443278019797 0.3564416360459859The reason for the non-consistency in your code is this:
On the line test_stat, test_p = stats.ttest_ind(men, women) you accepted the default setting that the t-test is to be computed by the equal-variances assumption. So the computation that scipy.stats gives you is a pure equal-variance t-test. That is described in the documentation for scipy.stats.ttest_ind
In your own code, you followed the Welch test in general: you computed the estimate of the mean and its standard error for the men and women separately, and computed the t-statistic in that way.
You did deviate from the Welch test in one place: the degrees-of-freedom computation. The degrees of freedom should be approximated with the formula I entered in the code (and linked to above), but you used the computation applicable under equal-variance assumptions.
IF you want more details on how to compute these statistics, or why they are defined as they are, or why your code is not what you expected it, I suggest you check out https://stats.stackexchange.com/ and https://datascience.stackexchange.com/ that are more appropriate for statistics questions, in comparison to https://stackoverflow.com/ that is more about the programming. Both those communities are proficient in python as well, so they should be able to help you out perfectly good.