[SciPy-Dev] Independent T-tests with unequal variances

Tue May 1 03:08:50 EDT 2012

Hi Gavin,

This is interesting. Thanks for sharing. Do you think that you could
submit this as a pull request on Github? A diff on a mailing list has
much more chances of getting lost in the flow.

Thanks a lot,

Gael

On Mon, Apr 30, 2012 at 06:11:47PM -0700, Junkshops wrote:
> Well, this is embarrassing. I appear to have sent the wrong diff
> with a sign error in it. The correct one is attached.

> My apologies,

> Gavin

> On 4/30/2012 5:12 PM, Junkshops wrote:
> >Hello all,

> >I hope, as an utter newb poking my nose into this list, that I'm
> >not giving the author of Miss Manner's Book of Netiquette the
> >vapours.

> >This is a follow up to Deniz Turgut's recent email:
> >http://article.gmane.org/gmane.comp.python.scientific.devel/16291/

> >"There is also a form of t-test for independent samples with
> >different variances, also known as Welch's t-test [2]. I think it
> >is better to include the 'identical variance' assumption in the
> >doc to avoid confusion."

> >I was just caught by this problem and heartily agree with Deniz's
> >views. However, I didn't see any explicit plans to add Welch's
> >test to scipy/stats/stats.py, and I needed an implementation of
> >the test and so implemented it. A diff against scipy 0.10.1 is
> >attached if anyone might find it useful.

> >Cheers,

> >Gavin

> 214c214
> <            'ttest_1samp', 'ttest_ind', 'ttest_rel',
> ---
> >            'ttest_1samp', 'ttest_ind', 'ttest_ind_uneq_var', 'ttest_rel',
> 2901c2901
> <     """Calculates the T-test for the means of TWO INDEPENDENT samples of scores.
> ---
> >     """Calculates the T-test for the means of TWO INDEPENDENT samples of scores with equal variances.
> 2903c2903
> <     This is a two-sided test for the null hypothesis that 2 independent samples
> ---
> >     This is a two-sided test for the null hypothesis that 2 independent samples with equal variances
> 2926c2926
> <     We can use this test, if we observe two independent samples from
> ---
> >     We can use this test, if we observe two independent samples with equal variances from
> 2972c2972
> < 
> ---

> 2988c2988
> < 
> ---

> 2990a2991,3116
> > def ttest_ind_uneq_var(a, b, axis=0):
> >     """Calculates the T-test for the means of TWO INDEPENDENT samples of scores with 
> >     unequal variances (or situations where it is unknown if the variances are equal).

> >     This is a two-sided test for the null hypothesis that 2 independent samples with unequal variances
> >     have identical average (expected) values.

> >     Parameters
> >     ----------
> >     a, b : sequence of ndarrays
> >         The arrays must have the same shape, except in the dimension
> >         corresponding to `axis` (the first, by default).
> >     axis : int, optional
> >         Axis can equal None (ravel array first), or an integer (the axis
> >         over which to operate on a and b).

> >     Returns
> >     -------
> >     t : float or array
> >         t-statistic
> >     prob : float or array
> >         two-tailed p-value

> >     Notes
> >     -----

> >     We can use this test, if we observe two independent samples from
> >     the same or different population, e.g. exam scores of boys and
> >     girls or of two ethnic groups. The test measures whether the
> >     average (expected) value differs significantly across samples. If
> >     we observe a large p-value, for example larger than 0.05 or 0.1,
> >     then we cannot reject the null hypothesis of identical average scores.
> >     If the p-value is smaller than the threshold, e.g. 1%, 5% or 10%,
> >     then we reject the null hypothesis of equal averages.

> >     References
> >     ----------

> >        http://en.wikipedia.org/wiki/T-test#Independent_two-sample_t-test

> >     Examples
> >     --------

> >     >>> from scipy import stats
> >     >>> import numpy as np

> >     >>> #fix seed to get the same result
> >     >>> np.random.seed(12345678)

> >     test with sample with identical means

> >     >>> rvs1 = stats.norm.rvs(loc=5,scale=10,size=500)
> >     >>> rvs2 = stats.norm.rvs(loc=5,scale=10,size=500)
> >     >>> stats.ttest_ind(rvs1,rvs2)
> >     (0.26833823296239279, 0.78849443369564765)
> >     >>> stats.ttest_ind_uneq_var(rvs1, rvs2)
> >     (0.26833823296239279, 0.78849419539158605)

> >     ttest_ind underestimates p for unequal variances

> >     >>> rvs3 = stats.norm.rvs(loc=5, scale=20, size=500)
> >     >>> stats.ttest_ind(rvs1, rvs3)
> >     (-0.46580283298287162, 0.64145827413436174)
> >     >>> stats.ttest_ind_uneq_var(rvs1, rvs3)
> >     (-0.46580283298287162, 0.64149552307593671)

> >     >>> rvs4 = stats.norm.rvs(loc=5, scale=20, size=100)
> >     >>> stats.ttest_ind(rvs1, rvs4)
> >     (-0.99882539442782481, 0.31828327091038955)
> >     >>> stats.ttest_ind_uneq_var(rvs1, rvs4)
> >     (-0.69712570584654099, 0.48711638692035597)

> >     test with sample with different means

> >     >>> rvs5 = stats.norm.rvs(loc=8, scale=10, size=500)
> >     >>> stats.ttest_ind(rvs1, rvs5)
> >     (-4.130511725493573, 3.922607411074624e-05)
> >     >>> stats.ttest_ind_uneq_var(rvs1, rvs5)
> >     (-4.130511725493573, 3.9209626240360421e-05)

> >     >>> rvs6 = stats.norm.rvs(loc=8, scale=20, size=500)
> >     >>> stats.ttest_ind(rvs1, rvs6)
> >     (-3.8383088416156559, 0.00013167799566923922)
> >     >>> stats.ttest_ind_uneq_var(rvs1, rvs6)
> >     (-3.8383088416156559, 0.00013475714831652827)

> >     >>> rvs7 = stats.norm.rvs(loc=8, scale=20, size=100)
> >     >>> stats.ttest_ind(rvs1, rvs7)
> >     (-0.79821473077740479, 0.42506275883963907)
> >     >>> stats.ttest_ind_uneq_var(rvs1, rvs7)
> >     (-0.51902756162811092, 0.60475596772293294)

> >     """

> >     a, b, axis = _chk2_asarray(a, b, axis)

> >     v1 = np.var(a, axis, ddof = 1)
> >     v2 = np.var(b, axis, ddof = 1)
> >     n1 = a.shape[axis]
> >     n2 = b.shape[axis]
> >     vn1 = v1 / n1
> >     vn2 = v2 / n2
> >     df = ((vn1 + vn2)**2) / ((vn1**2) / (n1 + 1) + (vn2**2) / (n2 + 1)) - 2

> >     d = np.mean(a, axis) - np.mean(b, axis)

> >     t = d / np.sqrt(vn1 + vn2)
> >     # why is this defining t = 1 when the means and variances are equal...?
> >     #t = np.where((d==0)*(svar==0), 1.0, t) #define t=0/0 = 0, identical means 
> >     #don't think this is really necessary, but returns t as array if not
> >     t = np.where((d==0), 0.0, t) #define t=0/0 = 0, identical means

> >     prob = distributions.t.sf(np.abs(t), df) * 2 #use np.abs to get upper tail

> >     #distributions.t.sf currently does not propagate nans
> >     #this can be dropped, if distributions.t.sf propagates nans
> >     #if this is removed, then prob = prob[()] needs to be removed
> >     prob = np.where(np.isnan(t), np.nan, prob)

> >     if t.ndim == 0:
> >         t = t[()]
> >         prob = prob[()]

> >     return t, prob

> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev

-- 
    Gael Varoquaux
    Researcher, INRIA Parietal
    Laboratoire de Neuro-Imagerie Assistee par Ordinateur
    NeuroSpin/CEA Saclay , Bat 145, 91191 Gif-sur-Yvette France
    Phone:  ++ 33-1-69-08-79-68
    http://gael-varoquaux.info            http://twitter.com/GaelVaroquaux