CODEWITHSUNDEEP
pythonnumpystatisticsscipy
Keith
Keith

Reputation: 4914

scipy.stats.ttest_ind without array (python)

I have done a number of calculations to estimate μ, σ and N for my two samples. Due to a number of approximations I don't have the arrays that are expected as input to scipy.stats.ttest_ind. Unless I am mistaken I only need μ, σ and N to do a welch's t test. Is there a way to do this in python?

Upvotes: 3

Views: 3625

Answers (3)

Josef
Josef

Reputation: 22897

As an update

The function is now available in scipy.stats, since version 0.16.0

http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.stats.ttest_ind_from_stats.html

scipy.stats.ttest_ind_from_stats(mean1, std1, nobs1, mean2, std2, nobs2, equal_var=True)
T-test for means of two independent samples from descriptive statistics.

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

Upvotes: 2

rroowwllaanndd
rroowwllaanndd

Reputation: 3958

Here’s a straightforward implementation based on this: 

import scipy.stats as stats
import numpy as np

def welch_t_test(mu1, s1, N1, mu2, s2, N2):
  # Construct arrays to make calculations more succint.
  N_i = np.array([N1, N2])
  dof_i = N_i - 1
  v_i = np.array([s1, s2]) ** 2
  # Calculate t-stat, degrees of freedom, use scipy to find p-value.
  t = (mu1 - mu2) / np.sqrt(np.sum(v_i / N_i))
  dof = (np.sum(v_i / N_i) ** 2) / np.sum((v_i ** 2) / ((N_i ** 2) * dof_i))
  p = stats.distributions.t.sf(np.abs(t), dof) * 2
  return t, p

It yields virtually identical results:

sample1 = np.random.rand(10)
sample2 = np.random.rand(15)
result_test = welch_t_test(np.mean(sample1), np.std(sample1, ddof=1), sample1.size,
                           np.mean(sample2), np.std(sample2, ddof=1), sample2.size)
result_scipy = stats.ttest_ind(sample1, sample2,equal_var=False)
np.allclose(result_test, result_scipy)
True

Upvotes: 3

Josef
Josef

Reputation: 22897

I have written t-test and z-test functions that take the summary statistics for statsmodels.

Those were intended mainly as internal shortcuts to avoid code duplication, and are not well documented.

For example http://statsmodels.sourceforge.net/devel/generated/statsmodels.stats.weightstats._tstat_generic.html

The list of related functions is here: http://statsmodels.sourceforge.net/devel/stats.html#basic-statistics-and-t-tests-with-frequency-weights

edit: in reply to comment

The function just does the core calculation, the actual calculation of the standard deviation of the difference under different assumptions is added to the calling method. https://github.com/statsmodels/statsmodels/blob/master/statsmodels/stats/weightstats.py#L713

edit

Here is an example how to use the methods of the CompareMeans class that includes the t-test based on summary statistics. We need to create a class that holds the relevant summary statistic as attributes. At the end there is a function that just wraps the relevant calls.

"""
Created on Wed Jul 23 05:47:34 2014
Author: Josef Perktold
License: BSD-3

"""

import numpy as np
from scipy import stats
from statsmodels.stats.weightstats import CompareMeans, ttest_ind


class SummaryStats(object):

    def __init__(self, nobs, mean, std):
        self.nobs = nobs
        self.mean = mean
        self.std = std
        self._var = std**2

np.random.seed(123)
nobs = 20
x1 = 1 + np.random.randn(nobs)
x2 = 1 + 1.5 * np.random.randn(nobs)

print stats.ttest_ind(x1, x2, equal_var=False)
print ttest_ind(x1, x2, usevar='unequal')

s1 = SummaryStats(x1.shape[0], x1.mean(0), x1.std(0))
s2 = SummaryStats(x2.shape[0], x2.mean(0), x2.std(0))

print CompareMeans(s1, s2).ttest_ind(usevar='unequal')



def ttest_ind_summ(summ1, summ2, usevar='unequal'):
    """t-test for equality of means based on summary statistic

    Parameters
    ----------
    summ1, summ2 : tuples of (nobs, mean, std)
        summary statistic for the two samples

    """

    s1 = SummaryStats(*summ1)
    s2 = SummaryStats(*summ2)
    return CompareMeans(s1, s2).ttest_ind(usevar=usevar)

print ttest_ind_summ((x1.shape[0], x1.mean(0), x1.std(0)),
                     (x2.shape[0], x2.mean(0), x2.std(0)), 
                     usevar='unequal')

''' result
(array(1.1590347327654558), 0.25416326823881513)
(1.1590347327654555, 0.25416326823881513, 35.573591346616553)
(1.1590347327654558, 0.25416326823881513, 35.57359134661656)
(1.1590347327654558, 0.25416326823881513, 35.57359134661656)
'''

Upvotes: 1

Related Questions