Last active
February 2, 2021 10:00
-
-
Save jdmonaco/5922991 to your computer and use it in GitHub Desktop.
Welch's t-test for two samples, not assuming equal sample size or variance. Requires Python, NumPy, and SciPy.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from collections import namedtuple | |
import numpy as np | |
import scipy.stats as st | |
TtestResults = namedtuple("Ttest", "T p") | |
def t_welch(x, y, tails=2): | |
""" | |
Welch's t-test for two unequal-size samples, not assuming equal variances | |
""" | |
assert tails in (1,2), "invalid: tails must be 1 or 2, found %s"%str(tails) | |
x, y = np.asarray(x), np.asarray(y) | |
nx, ny = x.size, y.size | |
vx, vy = x.var(ddof=1), y.var(ddof=1) | |
df = ((vx/nx + vy/ny)**2 / # Welch-Satterthwaite equation | |
((vx/nx)**2 / (nx - 1) + (vy/ny)**2 / (ny - 1))) | |
t_obs = (x.mean() - y.mean()) / np.sqrt(vx/nx + vy/ny) | |
p_value = tails * st.t.sf(abs(t_obs), df) | |
return TtestResults(t_obs, p_value) |
Thanks @jonkrohn! It's annoying because I pulled this from a larger stats module I wrote that had ddof=1
in every other std/var call. I somehow neglected these instances. And I was apparently under the impression in 2013 that dof's could only be integers, which is obviously not right (but at least the int
cast produced a conservative bias).
Anyway, I've updated this gist (and the corresponding module, which I do still use). Thanks again.
Nice, looks perfect now :)
I checked out your website, btw -- very cool! I did a neuroscience PhD and so naturally get jealous when I see fascinating work like yours by people who stayed on the academic course.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hi @jdmonaco! I was using this function of yours and noticed that the results are very close to being spot on.
After doing some digging around, I believe the sample-specific variance calculations should be adjusted slightly to account for a single degree of freedom because NumPy's
var()
method defaults to zero. That is:vx, vy = x.var(ddof=1), y.var(ddof=1)
Further, the degrees of freedom for calculating the p-value could optionally be executed with greater precision by dropping the
int()
in the Welch-Satterthwaite equation.