I am comparing stats.ttest_ind() vs “manual” computation of the same test, and get different results.
import numpy as np import pandas as pd import scipy.stats as stats import math
stats.ttest_ind() method:
#generate data
np.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 groups
men = df.age[df.sex == 'male']
women = df.age[df.sex == 'female']
#run t-test
test_stat, test_p = stats.ttest_ind(men, women)
print(test_stat, test_p)
Out:
-0.9265613940505325 0.355282312357339
Manual method:
#mean men_mean, women_mean = men.mean(), women.mean() #standard deviation men_sd, women_sd = men.std(ddof=1), women.std(ddof=1) #standard error men_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 women se_diff = math.sqrt(men_se**2.0 + women_se**2.0) #t-stat t_stat = (men_mean - women_mean) / se_diff #degrees of freedom df = men_n + women_n - 2 #critical value alpha = 0.05 cv = stats.t.ppf(1.0 - alpha, df) # p-value p = (1 - stats.t.cdf(abs(t_stat), df)) * 2 print(t_stat, cv, p)
Out:
-0.9244538916746341 0.3563753194455255
We 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 np
import pandas as pd
import scipy.stats as stats
import math
#generate data
np.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 groups
men = df.age[df.sex == 'male']
women = df.age[df.sex == 'female']
#run t-test
############################### CHANGED THE ROW BELOW HERE
test_stat, test_p = stats.ttest_ind(men, women,equal_var=False)
print(test_stat, test_p)
#mean
men_mean, women_mean = men.mean(), women.mean()
#standard deviation
men_sd, women_sd = men.std(ddof=1), women.std(ddof=1)
#standard error
men_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 women
se_diff = math.sqrt(men_se**2.0 + women_se**2.0)
#t-stat
t_stat = (men_mean - women_mean) / se_diff
#degrees of freedom
############################### CHANGED THE ROW BELOW HERE
df = (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 value
alpha = 0.05
cv = stats.t.ppf(1.0 - alpha, df)
# p-value
p = (1 - stats.t.cdf(abs(t_stat), df)) * 2
print(t_stat, cv, p)
and it outputs
-0.9244538916746341 0.356441636045986 -0.9244538916746341 1.6530443278019797 0.3564416360459859
The 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.