# Tracking down the assumptions made by SciPy's ttest_ind() function

I'm trying to write my own Python code to compute t-statistics and p-values for one and two tailed independent t tests. I can use the normal approximation, but for the moment I am trying to just use the t-distribution. I've been unsuccessful in matching the results of SciPy's stats library on my test data. I could use a fresh pair of eyes to see if I'm just making a dumb mistake somewhere.

Note, this isn't so much of a coding question as it is a "why isn't this computation yielding the right t-stat?" I give the code for completeness, but don't expect any software advice. Just help understanding why this isn't right.

My code:

import numpy as np
import scipy.stats as st

def compute_t_stat(pop1,pop2):

num1 = pop1.shape[0]; num2 = pop2.shape[0];

# The formula for t-stat when population variances differ.
t_stat = (np.mean(pop1) - np.mean(pop2))/np.sqrt( np.var(pop1)/num1 + np.var(pop2)/num2 )

# ADDED: The Welch-Satterthwaite degrees of freedom.
df = ((np.var(pop1)/num1 + np.var(pop2)/num2)**(2.0))/(   (np.var(pop1)/num1)**(2.0)/(num1-1) +  (np.var(pop2)/num2)**(2.0)/(num2-1) )

# Am I computing this wrong?
# It should just come from the CDF like this, right?
# The extra parameter is the degrees of freedom.

one_tailed_p_value = 1.0 - st.t.cdf(t_stat,df)
two_tailed_p_value = 1.0 - ( st.t.cdf(np.abs(t_stat),df) - st.t.cdf(-np.abs(t_stat),df) )

# Computing with SciPy's built-ins
# My results don't match theirs.
t_ind, p_ind = st.ttest_ind(pop1, pop2)

return t_stat, one_tailed_p_value, two_tailed_p_value, t_ind, p_ind


Update:

After reading a bit more on the Welch's t-test, I saw that I should be using the Welch-Satterthwaite formula to calculate degrees of freedom. I updated the code above to reflect this.

With the new degrees of freedom, I get a closer result. My two-sided p-value is off by about 0.008 from the SciPy version's... but this is still much too big an error so I must still be doing something incorrect (or SciPy distribution functions are very bad, but it's hard to believe they are only accurate to 2 decimal places).

Second update:

While continuing to try things, I thought maybe SciPy's version automatically computes the Normal approximation to the t-distribution when the degrees of freedom are high enough (roughly > 30). So I re-ran my code using the Normal distribution instead, and the computed results are actually further away from SciPy's than when I use the t-distribution.

Bonus question :)

Also, the t-statistic is negative. I was just wondering what this means for the one-sided t-test. Does this typically mean that I should be looking in the negative axis direction for the test? In my test data, population 1 is a control group who did not receive a certain employment training program. Population 2 did receive it, and the measured data are wage differences before/after treatment.

So I have some reason to think that the mean for population 2 will be larger. But from a statistical theory point of view, it doesn't seem right to concoct a test this way. How could I have known to check (for the one-sided test) in the negative direction without relying on subjective knowledge about the data? Or is this just one of those frequentist things that, while not philosophically rigorous, needs to be done in practice?

-
Maybe SciPy calculates Welch's t-test -- SciPy's documentation doesn't specify... – Cyan Apr 5 '12 at 21:57
The formula I am using in my calculation is the same as the Welch's t-statistic. To my knowledge, this is the "standard" thing to do when the sample sizes and population variances are allowed to be different, correct? – EMS Apr 5 '12 at 21:59
Don't you need to take the square of the (current) numerator in the calculation of the degrees of freedom? Also, with virtually no code changes, there are much safer ways of calculating the $p$-values. The way it's currently implemented is extremely susceptible to massive error due to cancellation. – cardinal Apr 6 '12 at 3:26
(1) Check the documentation of numpy.var. The version I saw seems to indicate that the MLE estimate is calculated by default instead of the unbiased estimate. To get the unbiased estimate one needs to call it with the optional ddof=1. (2) For the upper-tail $p$-value, use the symmetry of the $t$-distribution, i.e., one_tailed_p_value = st.t.cdf(-t_stat,df) and (3) for the two-tailed $p$-value, do something similar: two_tailed_p_value = 2*st.t.cdf(-np.abs(t_stat),df). – cardinal Apr 6 '12 at 3:44
I don't think of it as all that trivial, in the sense that there is often a sizable gap between having a mathematical formula for something at hand and knowing a safe and efficient way of computing it. It's one of those things where it's nice to have a large body of knowledge already available, because it would take a virtual eternity to learn such tricks, one-by-one, all on your own. :) – cardinal Apr 6 '12 at 4:15