3

What is the most straightforward way to perform a t-test in python and to include CIs of the difference? I've seen various posts but everything is different and when I tried to calculate the CIs myself it seemed slightly wrong... Here:

import numpy as np
from scipy import stats

g1 = np.array([48.7107107,
36.8587287,
67.7129929,
39.5538852,
35.8622661])
g2 = np.array([62.4993857,
49.7434833,
67.7516511,
54.3585559,
71.0933957])

m1, m2 = np.mean(g1), np.mean(g2)
dof = (len(g1)-1) + (len(g2)-1)

MSE = (np.var(g1) + np.var(g2)) / 2

stderr_diffs = np.sqrt((2 * MSE)/len(g1))

tcl = stats.t.ppf([.975], dof)

lower_limit = (m1-m2) - (tcl) * (stderr_diffs)
upper_limit = (m1-m2) + (tcl) * (stderr_diffs)

print(lower_limit, upper_limit)

returns:

[-30.12845447] [-0.57070077]

However, when I run the same test in SPSS, although I have the same t and p values, the CIs are -31.87286, 1.17371, and this is also the case in R. I can't seem to find the correct way to do this and would appreciate some help.

fffrost
  • 1,659
  • 1
  • 21
  • 36

1 Answers1

3

You're subtracting 1 when you compute dof, but when you compute the variance you're not using the sample variance:

MSE = (np.var(g1) + np.var(g2)) / 2

should be

MSE = (np.var(g1, ddof=1) + np.var(g2, ddof=1)) / 2

which gives me

[-31.87286426] [ 1.17370902]

That said, instead of doing the manual implementation, I'd probably use statsmodels' CompareMeans:

In [105]: import statsmodels.stats.api as sms

In [106]: r = sms.CompareMeans(sms.DescrStatsW(g1), sms.DescrStatsW(g2))

In [107]: r.tconfint_diff()
Out[107]: (-31.872864255548553, 1.1737090155485568)

(really we should be using a DataFrame here, not an ndarray, but I'm lazy).

Remember though that you're going to want to consider what assumption you want to make about the variance:

In [110]: r.tconfint_diff(usevar='pooled')
Out[110]: (-31.872864255548553, 1.1737090155485568)

In [111]: r.tconfint_diff(usevar='unequal')
Out[111]: (-32.28794665832114, 1.5887914183211436)

and if your g1 and g2 are representative, the assumption of equal variance might not be a good one.

DSM
  • 342,061
  • 65
  • 592
  • 494
  • Thanks for the answer. I get the following when I import statsmodels.stats.api as sms: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead. from pandas.core import datetools Is there anything I can do about this? – fffrost Feb 10 '18 at 16:32
  • Unfortunately, not really -- that's a problem on the statsmodels side, where they haven't updated yet. In the meantime you could add `import warnings` and `warnings.simplefilter(action='ignore', category=FutureWarning)` before the statsmodels import if you just want to silence it (note: untested, so I might have got the syntax wrong, but 'python suppress futurewarning' or something would probably bring the right way up.) – DSM Feb 11 '18 at 16:23