Source code for statsmodels.stats.weightstats
"""Tests and descriptive statistics with weights
Created on 2010-09-18
Author: josef-pktd
License: BSD (3-clause)
References
----------
SPSS manual
SAS manual
This follows in large parts the SPSS manual, which is largely the same as
the SAS manual with different, simpler notation.
Freq, Weight in SAS seems redundant since they always show up as product, SPSS
has only weights.
Notes
-----
This has potential problems with ddof, I started to follow numpy with ddof=0
by default and users can change it, but this might still mess up the t-tests,
since the estimates for the standard deviation will be based on the ddof that
the user chooses.
- fixed ddof for the meandiff ttest, now matches scipy.stats.ttest_ind
Note: scipy has now a separate, pooled variance option in ttest, but I have not
compared yet.
"""
import numpy as np
from scipy import stats
from statsmodels.tools.decorators import cache_readonly
[docs]
class DescrStatsW:
"""
Descriptive statistics and tests with weights for case weights
Assumes that the data is 1d or 2d with (nobs, nvars) observations in rows,
variables in columns, and that the same weight applies to each column.
If degrees of freedom correction is used, then weights should add up to the
number of observations. ttest also assumes that the sum of weights
corresponds to the sample size.
This is essentially the same as replicating each observations by its
weight, if the weights are integers, often called case or frequency weights.
Parameters
----------
data : array_like, 1-D or 2-D
dataset
weights : None or 1-D ndarray
weights for each observation, with same length as zero axis of data
ddof : int
default ddof=0, degrees of freedom correction used for second moments,
var, std, cov, corrcoef.
However, statistical tests are independent of `ddof`, based on the
standard formulas.
Examples
--------
>>> import numpy as np
>>> np.random.seed(0)
>>> x1_2d = 1.0 + np.random.randn(20, 3)
>>> w1 = np.random.randint(1, 4, 20)
>>> d1 = DescrStatsW(x1_2d, weights=w1)
>>> d1.mean
array([ 1.42739844, 1.23174284, 1.083753 ])
>>> d1.var
array([ 0.94855633, 0.52074626, 1.12309325])
>>> d1.std_mean
array([ 0.14682676, 0.10878944, 0.15976497])
>>> tstat, pval, df = d1.ttest_mean(0)
>>> tstat; pval; df
array([ 9.72165021, 11.32226471, 6.78342055])
array([ 1.58414212e-12, 1.26536887e-14, 2.37623126e-08])
44.0
>>> tstat, pval, df = d1.ttest_mean([0, 1, 1])
>>> tstat; pval; df
array([ 9.72165021, 2.13019609, 0.52422632])
array([ 1.58414212e-12, 3.87842808e-02, 6.02752170e-01])
44.0
# if weights are integers, then asrepeats can be used
>>> x1r = d1.asrepeats()
>>> x1r.shape
...
>>> stats.ttest_1samp(x1r, [0, 1, 1])
...
"""
def __init__(self, data, weights=None, ddof=0):
self.data = np.asarray(data)
if weights is None:
self.weights = np.ones(self.data.shape[0])
else:
self.weights = np.asarray(weights).astype(float)
# TODO: why squeeze?
if len(self.weights.shape) > 1 and len(self.weights) > 1:
self.weights = self.weights.squeeze()
self.ddof = ddof
@cache_readonly
def sum_weights(self):
"""Sum of weights"""
return self.weights.sum(0)
@cache_readonly
def nobs(self):
"""alias for number of observations/cases, equal to sum of weights
"""
return self.sum_weights
@cache_readonly
def sum(self):
"""weighted sum of data"""
return np.dot(self.data.T, self.weights)
@cache_readonly
def mean(self):
"""weighted mean of data"""
return self.sum / self.sum_weights
@cache_readonly
def demeaned(self):
"""data with weighted mean subtracted"""
return self.data - self.mean
@cache_readonly
def sumsquares(self):
"""weighted sum of squares of demeaned data"""
return np.dot((self.demeaned ** 2).T, self.weights)
# need memoize instead of cache decorator
[docs]
def var_ddof(self, ddof=0):
"""variance of data given ddof
Parameters
----------
ddof : int, float
degrees of freedom correction, independent of attribute ddof
Returns
-------
var : float, ndarray
variance with denominator ``sum_weights - ddof``
"""
return self.sumsquares / (self.sum_weights - ddof)
[docs]
def std_ddof(self, ddof=0):
"""standard deviation of data with given ddof
Parameters
----------
ddof : int, float
degrees of freedom correction, independent of attribute ddof
Returns
-------
std : float, ndarray
standard deviation with denominator ``sum_weights - ddof``
"""
return np.sqrt(self.var_ddof(ddof=ddof))
@cache_readonly
def var(self):
"""variance with default degrees of freedom correction
"""
return self.sumsquares / (self.sum_weights - self.ddof)
@cache_readonly
def _var(self):
"""variance without degrees of freedom correction
used for statistical tests with controlled ddof
"""
return self.sumsquares / self.sum_weights
@cache_readonly
def std(self):
"""standard deviation with default degrees of freedom correction
"""
return np.sqrt(self.var)
@cache_readonly
def cov(self):
"""weighted covariance of data if data is 2 dimensional
assumes variables in columns and observations in rows
uses default ddof
"""
cov_ = np.dot(self.weights * self.demeaned.T, self.demeaned)
cov_ /= self.sum_weights - self.ddof
return cov_
@cache_readonly
def corrcoef(self):
"""weighted correlation with default ddof
assumes variables in columns and observations in rows
"""
return self.cov / self.std / self.std[:, None]
@cache_readonly
def std_mean(self):
"""standard deviation of weighted mean
"""
std = self.std
if self.ddof != 0:
# ddof correction, (need copy of std)
std = std * np.sqrt(
(self.sum_weights - self.ddof) / self.sum_weights
)
return std / np.sqrt(self.sum_weights - 1)
[docs]
def quantile(self, probs, return_pandas=True):
"""
Compute quantiles for a weighted sample.
Parameters
----------
probs : array_like
A vector of probability points at which to calculate the
quantiles. Each element of `probs` should fall in [0, 1].
return_pandas : bool
If True, return value is a Pandas DataFrame or Series.
Otherwise returns a ndarray.
Returns
-------
quantiles : Series, DataFrame, or ndarray
If `return_pandas` = True, returns one of the following:
* data are 1d, `return_pandas` = True: a Series indexed by
the probability points.
* data are 2d, `return_pandas` = True: a DataFrame with
the probability points as row index and the variables
as column index.
If `return_pandas` = False, returns an ndarray containing the
same values as the Series/DataFrame.
Notes
-----
To compute the quantiles, first, the weights are summed over
exact ties yielding distinct data values y_1 < y_2 < ..., and
corresponding weights w_1, w_2, .... Let s_j denote the sum
of the first j weights, and let W denote the sum of all the
weights. For a probability point p, if pW falls strictly
between s_j and s_{j+1} then the estimated quantile is
y_{j+1}. If pW = s_j then the estimated quantile is (y_j +
y_{j+1})/2. If pW < p_1 then the estimated quantile is y_1.
References
----------
SAS documentation for weighted quantiles:
https://support.sas.com/documentation/cdl/en/procstat/63104/HTML/default/viewer.htm#procstat_univariate_sect028.htm
"""
import pandas as pd
probs = np.asarray(probs)
probs = np.atleast_1d(probs)
if self.data.ndim == 1:
rslt = self._quantile(self.data, probs)
if return_pandas:
rslt = pd.Series(rslt, index=probs)
else:
rslt = []
for vec in self.data.T:
rslt.append(self._quantile(vec, probs))
rslt = np.column_stack(rslt)
if return_pandas:
columns = ["col%d" % (j + 1) for j in range(rslt.shape[1])]
rslt = pd.DataFrame(data=rslt, columns=columns, index=probs)
if return_pandas:
rslt.index.name = "p"
return rslt
def _quantile(self, vec, probs):
# Helper function to calculate weighted quantiles for one column.
# Follows definition from SAS documentation.
# Returns ndarray
import pandas as pd
# Aggregate over ties
df = pd.DataFrame(index=np.arange(len(self.weights)))
df["weights"] = self.weights
df["vec"] = vec
dfg = df.groupby("vec").agg("sum")
weights = dfg.values[:, 0]
values = np.asarray(dfg.index)
cweights = np.cumsum(weights)
totwt = cweights[-1]
targets = probs * totwt
ii = np.searchsorted(cweights, targets)
rslt = values[ii]
# Exact hits
jj = np.flatnonzero(np.abs(targets - cweights[ii]) < 1e-10)
jj = jj[ii[jj] < len(cweights) - 1]
rslt[jj] = (values[ii[jj]] + values[ii[jj] + 1]) / 2
return rslt
[docs]
def tconfint_mean(self, alpha=0.05, alternative="two-sided"):
"""two-sided confidence interval for weighted mean of data
If the data is 2d, then these are separate confidence intervals
for each column.
Parameters
----------
alpha : float
significance level for the confidence interval, coverage is
``1-alpha``
alternative : str
This specifies the alternative hypothesis for the test that
corresponds to the confidence interval.
The alternative hypothesis, H1, has to be one of the following
'two-sided': H1: mean not equal to value (default)
'larger' : H1: mean larger than value
'smaller' : H1: mean smaller than value
Returns
-------
lower, upper : floats or ndarrays
lower and upper bound of confidence interval
Notes
-----
In a previous version, statsmodels 0.4, alpha was the confidence
level, e.g. 0.95
"""
# TODO: add asymmetric
dof = self.sum_weights - 1
ci = _tconfint_generic(
self.mean, self.std_mean, dof, alpha, alternative
)
return ci
[docs]
def zconfint_mean(self, alpha=0.05, alternative="two-sided"):
"""two-sided confidence interval for weighted mean of data
Confidence interval is based on normal distribution.
If the data is 2d, then these are separate confidence intervals
for each column.
Parameters
----------
alpha : float
significance level for the confidence interval, coverage is
``1-alpha``
alternative : str
This specifies the alternative hypothesis for the test that
corresponds to the confidence interval.
The alternative hypothesis, H1, has to be one of the following
'two-sided': H1: mean not equal to value (default)
'larger' : H1: mean larger than value
'smaller' : H1: mean smaller than value
Returns
-------
lower, upper : floats or ndarrays
lower and upper bound of confidence interval
Notes
-----
In a previous version, statsmodels 0.4, alpha was the confidence
level, e.g. 0.95
"""
return _zconfint_generic(self.mean, self.std_mean, alpha, alternative)
[docs]
def ttest_mean(self, value=0, alternative="two-sided"):
"""ttest of Null hypothesis that mean is equal to value.
The alternative hypothesis H1 is defined by the following
- 'two-sided': H1: mean not equal to value
- 'larger' : H1: mean larger than value
- 'smaller' : H1: mean smaller than value
Parameters
----------
value : float or array
the hypothesized value for the mean
alternative : str
The alternative hypothesis, H1, has to be one of the following:
- 'two-sided': H1: mean not equal to value (default)
- 'larger' : H1: mean larger than value
- 'smaller' : H1: mean smaller than value
Returns
-------
tstat : float
test statistic
pvalue : float
pvalue of the t-test
df : int or float
"""
# TODO: check direction with R, smaller=less, larger=greater
tstat = (self.mean - value) / self.std_mean
dof = self.sum_weights - 1
# TODO: use outsourced
if alternative == "two-sided":
pvalue = stats.t.sf(np.abs(tstat), dof) * 2
elif alternative == "larger":
pvalue = stats.t.sf(tstat, dof)
elif alternative == "smaller":
pvalue = stats.t.cdf(tstat, dof)
else:
raise ValueError("alternative not recognized")
return tstat, pvalue, dof
[docs]
def ttost_mean(self, low, upp):
"""test of (non-)equivalence of one sample
TOST: two one-sided t tests
null hypothesis: m < low or m > upp
alternative hypothesis: low < m < upp
where m is the expected value of the sample (mean of the population).
If the pvalue is smaller than a threshold, say 0.05, then we reject the
hypothesis that the expected value of the sample (mean of the
population) is outside of the interval given by thresholds low and upp.
Parameters
----------
low, upp : float
equivalence interval low < mean < upp
Returns
-------
pvalue : float
pvalue of the non-equivalence test
t1, pv1, df1 : tuple
test statistic, pvalue and degrees of freedom for lower threshold
test
t2, pv2, df2 : tuple
test statistic, pvalue and degrees of freedom for upper threshold
test
"""
t1, pv1, df1 = self.ttest_mean(low, alternative="larger")
t2, pv2, df2 = self.ttest_mean(upp, alternative="smaller")
return np.maximum(pv1, pv2), (t1, pv1, df1), (t2, pv2, df2)
[docs]
def ztest_mean(self, value=0, alternative="two-sided"):
"""z-test of Null hypothesis that mean is equal to value.
The alternative hypothesis H1 is defined by the following
'two-sided': H1: mean not equal to value
'larger' : H1: mean larger than value
'smaller' : H1: mean smaller than value
Parameters
----------
value : float or array
the hypothesized value for the mean
alternative : str
The alternative hypothesis, H1, has to be one of the following
'two-sided': H1: mean not equal to value (default)
'larger' : H1: mean larger than value
'smaller' : H1: mean smaller than value
Returns
-------
tstat : float
test statistic
pvalue : float
pvalue of the z-test
Notes
-----
This uses the same degrees of freedom correction as the t-test in the
calculation of the standard error of the mean, i.e it uses
`(sum_weights - 1)` instead of `sum_weights` in the denominator.
See Examples below for the difference.
Examples
--------
z-test on a proportion, with 20 observations, 15 of those are our event
>>> import statsmodels.api as sm
>>> x1 = [0, 1]
>>> w1 = [5, 15]
>>> d1 = sm.stats.DescrStatsW(x1, w1)
>>> d1.ztest_mean(0.5)
(2.5166114784235836, 0.011848940928347452)
This differs from the proportions_ztest because of the degrees of
freedom correction:
>>> sm.stats.proportions_ztest(15, 20.0, value=0.5)
(2.5819888974716112, 0.009823274507519247).
We can replicate the results from ``proportions_ztest`` if we increase
the weights to have artificially one more observation:
>>> sm.stats.DescrStatsW(x1, np.array(w1)*21./20).ztest_mean(0.5)
(2.5819888974716116, 0.0098232745075192366)
"""
tstat = (self.mean - value) / self.std_mean
# TODO: use outsourced
if alternative == "two-sided":
pvalue = stats.norm.sf(np.abs(tstat)) * 2
elif alternative == "larger":
pvalue = stats.norm.sf(tstat)
elif alternative == "smaller":
pvalue = stats.norm.cdf(tstat)
return tstat, pvalue
[docs]
def ztost_mean(self, low, upp):
"""test of (non-)equivalence of one sample, based on z-test
TOST: two one-sided z-tests
null hypothesis: m < low or m > upp
alternative hypothesis: low < m < upp
where m is the expected value of the sample (mean of the population).
If the pvalue is smaller than a threshold, say 0.05, then we reject the
hypothesis that the expected value of the sample (mean of the
population) is outside of the interval given by thresholds low and upp.
Parameters
----------
low, upp : float
equivalence interval low < mean < upp
Returns
-------
pvalue : float
pvalue of the non-equivalence test
t1, pv1 : tuple
test statistic and p-value for lower threshold test
t2, pv2 : tuple
test statistic and p-value for upper threshold test
"""
t1, pv1 = self.ztest_mean(low, alternative="larger")
t2, pv2 = self.ztest_mean(upp, alternative="smaller")
return np.maximum(pv1, pv2), (t1, pv1), (t2, pv2)
[docs]
def get_compare(self, other, weights=None):
"""return an instance of CompareMeans with self and other
Parameters
----------
other : array_like or instance of DescrStatsW
If array_like then this creates an instance of DescrStatsW with
the given weights.
weights : None or array
weights are only used if other is not an instance of DescrStatsW
Returns
-------
cm : instance of CompareMeans
the instance has self attached as d1 and other as d2.
See Also
--------
CompareMeans
"""
if not isinstance(other, self.__class__):
d2 = DescrStatsW(other, weights)
else:
d2 = other
return CompareMeans(self, d2)
[docs]
def asrepeats(self):
"""get array that has repeats given by floor(weights)
observations with weight=0 are dropped
"""
w_int = np.floor(self.weights).astype(int)
return np.repeat(self.data, w_int, axis=0)
[docs]
def _tstat_generic(value1, value2, std_diff, dof, alternative, diff=0):
"""generic ttest based on summary statistic
The test statistic is :
tstat = (value1 - value2 - diff) / std_diff
and is assumed to be t-distributed with ``dof`` degrees of freedom.
Parameters
----------
value1 : float or ndarray
Value, for example mean, of the first sample.
value2 : float or ndarray
Value, for example mean, of the second sample.
std_diff : float or ndarray
Standard error of the difference value1 - value2
dof : int or float
Degrees of freedom
alternative : str
The alternative hypothesis, H1, has to be one of the following
* 'two-sided' : H1: ``value1 - value2 - diff`` not equal to 0.
* 'larger' : H1: ``value1 - value2 - diff > 0``
* 'smaller' : H1: ``value1 - value2 - diff < 0``
diff : float
value of difference ``value1 - value2`` under the null hypothesis
Returns
-------
tstat : float or ndarray
Test statistic.
pvalue : float or ndarray
P-value of the hypothesis test assuming that the test statistic is
t-distributed with ``df`` degrees of freedom.
"""
tstat = (value1 - value2 - diff) / std_diff
if alternative in ["two-sided", "2-sided", "2s"]:
pvalue = stats.t.sf(np.abs(tstat), dof) * 2
elif alternative in ["larger", "l"]:
pvalue = stats.t.sf(tstat, dof)
elif alternative in ["smaller", "s"]:
pvalue = stats.t.cdf(tstat, dof)
else:
raise ValueError("invalid alternative")
return tstat, pvalue
[docs]
def _tconfint_generic(mean, std_mean, dof, alpha, alternative):
"""generic t-confint based on summary statistic
Parameters
----------
mean : float or ndarray
Value, for example mean, of the first sample.
std_mean : float or ndarray
Standard error of the difference value1 - value2
dof : int or float
Degrees of freedom
alpha : float
Significance level for the confidence interval, coverage is
``1-alpha``.
alternative : str
The alternative hypothesis, H1, has to be one of the following
* 'two-sided' : H1: ``value1 - value2 - diff`` not equal to 0.
* 'larger' : H1: ``value1 - value2 - diff > 0``
* 'smaller' : H1: ``value1 - value2 - diff < 0``
Returns
-------
lower : float or ndarray
Lower confidence limit. This is -inf for the one-sided alternative
"smaller".
upper : float or ndarray
Upper confidence limit. This is inf for the one-sided alternative
"larger".
"""
if alternative in ["two-sided", "2-sided", "2s"]:
tcrit = stats.t.ppf(1 - alpha / 2.0, dof)
lower = mean - tcrit * std_mean
upper = mean + tcrit * std_mean
elif alternative in ["larger", "l"]:
tcrit = stats.t.ppf(alpha, dof)
lower = mean + tcrit * std_mean
upper = np.inf
elif alternative in ["smaller", "s"]:
tcrit = stats.t.ppf(1 - alpha, dof)
lower = -np.inf
upper = mean + tcrit * std_mean
else:
raise ValueError("invalid alternative")
return lower, upper
[docs]
def _zstat_generic(value1, value2, std_diff, alternative, diff=0):
"""generic (normal) z-test based on summary statistic
The test statistic is :
tstat = (value1 - value2 - diff) / std_diff
and is assumed to be normally distributed.
Parameters
----------
value1 : float or ndarray
Value, for example mean, of the first sample.
value2 : float or ndarray
Value, for example mean, of the second sample.
std_diff : float or ndarray
Standard error of the difference value1 - value2
alternative : str
The alternative hypothesis, H1, has to be one of the following
* 'two-sided' : H1: ``value1 - value2 - diff`` not equal to 0.
* 'larger' : H1: ``value1 - value2 - diff > 0``
* 'smaller' : H1: ``value1 - value2 - diff < 0``
diff : float
value of difference ``value1 - value2`` under the null hypothesis
Returns
-------
tstat : float or ndarray
Test statistic.
pvalue : float or ndarray
P-value of the hypothesis test assuming that the test statistic is
t-distributed with ``df`` degrees of freedom.
"""
zstat = (value1 - value2 - diff) / std_diff
if alternative in ["two-sided", "2-sided", "2s"]:
pvalue = stats.norm.sf(np.abs(zstat)) * 2
elif alternative in ["larger", "l"]:
pvalue = stats.norm.sf(zstat)
elif alternative in ["smaller", "s"]:
pvalue = stats.norm.cdf(zstat)
else:
raise ValueError("invalid alternative")
return zstat, pvalue
[docs]
def _zstat_generic2(value, std, alternative):
"""generic (normal) z-test based on summary statistic
The test statistic is :
zstat = value / std
and is assumed to be normally distributed with standard deviation ``std``.
Parameters
----------
value : float or ndarray
Value of a sample statistic, for example mean.
value2 : float or ndarray
Value, for example mean, of the second sample.
std : float or ndarray
Standard error of the sample statistic value.
alternative : str
The alternative hypothesis, H1, has to be one of the following
* 'two-sided' : H1: ``value1 - value2 - diff`` not equal to 0.
* 'larger' : H1: ``value1 - value2 - diff > 0``
* 'smaller' : H1: ``value1 - value2 - diff < 0``
Returns
-------
zstat : float or ndarray
Test statistic.
pvalue : float or ndarray
P-value of the hypothesis test assuming that the test statistic is
normally distributed.
"""
zstat = value / std
if alternative in ["two-sided", "2-sided", "2s"]:
pvalue = stats.norm.sf(np.abs(zstat)) * 2
elif alternative in ["larger", "l"]:
pvalue = stats.norm.sf(zstat)
elif alternative in ["smaller", "s"]:
pvalue = stats.norm.cdf(zstat)
else:
raise ValueError("invalid alternative")
return zstat, pvalue
[docs]
def _zconfint_generic(mean, std_mean, alpha, alternative):
"""generic normal-confint based on summary statistic
Parameters
----------
mean : float or ndarray
Value, for example mean, of the first sample.
std_mean : float or ndarray
Standard error of the difference value1 - value2
alpha : float
Significance level for the confidence interval, coverage is
``1-alpha``
alternative : str
The alternative hypothesis, H1, has to be one of the following
* 'two-sided' : H1: ``value1 - value2 - diff`` not equal to 0.
* 'larger' : H1: ``value1 - value2 - diff > 0``
* 'smaller' : H1: ``value1 - value2 - diff < 0``
Returns
-------
lower : float or ndarray
Lower confidence limit. This is -inf for the one-sided alternative
"smaller".
upper : float or ndarray
Upper confidence limit. This is inf for the one-sided alternative
"larger".
"""
if alternative in ["two-sided", "2-sided", "2s"]:
zcrit = stats.norm.ppf(1 - alpha / 2.0)
lower = mean - zcrit * std_mean
upper = mean + zcrit * std_mean
elif alternative in ["larger", "l"]:
zcrit = stats.norm.ppf(alpha)
lower = mean + zcrit * std_mean
upper = np.inf
elif alternative in ["smaller", "s"]:
zcrit = stats.norm.ppf(1 - alpha)
lower = -np.inf
upper = mean + zcrit * std_mean
else:
raise ValueError("invalid alternative")
return lower, upper
[docs]
class CompareMeans:
"""class for two sample comparison
The tests and the confidence interval work for multi-endpoint comparison:
If d1 and d2 have the same number of rows, then each column of the data
in d1 is compared with the corresponding column in d2.
Parameters
----------
d1, d2 : instances of DescrStatsW
Notes
-----
The result for the statistical tests and the confidence interval are
independent of the user specified ddof.
TODO: Extend to any number of groups or write a version that works in that
case, like in SAS and SPSS.
"""
def __init__(self, d1, d2):
"""assume d1, d2 hold the relevant attributes
"""
self.d1 = d1
self.d2 = d2
# assume nobs is available
# if not hasattr(self.d1, 'nobs'):
# d1.nobs1 = d1.sum_weights.astype(float) #float just to make sure
# self.nobs2 = d2.sum_weights.astype(float)
[docs]
@classmethod
def from_data(
cls, data1, data2, weights1=None, weights2=None, ddof1=0, ddof2=0
):
"""construct a CompareMeans object from data
Parameters
----------
data1, data2 : array_like, 1-D or 2-D
compared datasets
weights1, weights2 : None or 1-D ndarray
weights for each observation of data1 and data2 respectively,
with same length as zero axis of corresponding dataset.
ddof1, ddof2 : int
default ddof1=0, ddof2=0, degrees of freedom for data1,
data2 respectively.
Returns
-------
A CompareMeans instance.
"""
return cls(
DescrStatsW(data1, weights=weights1, ddof=ddof1),
DescrStatsW(data2, weights=weights2, ddof=ddof2),
)
[docs]
def summary(self, use_t=True, alpha=0.05, usevar="pooled", value=0):
"""summarize the results of the hypothesis test
Parameters
----------
use_t : bool, optional
if use_t is True, then t test results are returned
if use_t is False, then z test results are returned
alpha : float
significance level for the confidence interval, coverage is
``1-alpha``
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is
assumed to be the same. If ``unequal``, then the variance of
Welch ttest will be used, and the degrees of freedom are those
of Satterthwaite if ``use_t`` is True.
value : float
difference between the means under the Null hypothesis.
Returns
-------
smry : SimpleTable
"""
d1 = self.d1
d2 = self.d2
confint_percents = 100 - alpha * 100
if use_t:
tstat, pvalue, _ = self.ttest_ind(usevar=usevar, value=value)
lower, upper = self.tconfint_diff(alpha=alpha, usevar=usevar)
else:
tstat, pvalue = self.ztest_ind(usevar=usevar, value=value)
lower, upper = self.zconfint_diff(alpha=alpha, usevar=usevar)
if usevar == "pooled":
std_err = self.std_meandiff_pooledvar
else:
std_err = self.std_meandiff_separatevar
std_err = np.atleast_1d(std_err)
tstat = np.atleast_1d(tstat)
pvalue = np.atleast_1d(pvalue)
lower = np.atleast_1d(lower)
upper = np.atleast_1d(upper)
conf_int = np.column_stack((lower, upper))
params = np.atleast_1d(d1.mean - d2.mean - value)
title = "Test for equality of means"
yname = "y" # not used in params_frame
xname = ["subset #%d" % (ii + 1) for ii in range(tstat.shape[0])]
from statsmodels.iolib.summary import summary_params
return summary_params(
(None, params, std_err, tstat, pvalue, conf_int),
alpha=alpha,
use_t=use_t,
yname=yname,
xname=xname,
title=title,
)
@cache_readonly
def std_meandiff_separatevar(self):
# this uses ``_var`` to use ddof=0 for formula
d1 = self.d1
d2 = self.d2
return np.sqrt(d1._var / (d1.nobs - 1) + d2._var / (d2.nobs - 1))
@cache_readonly
def std_meandiff_pooledvar(self):
"""variance assuming equal variance in both data sets
"""
# this uses ``_var`` to use ddof=0 for formula
d1 = self.d1
d2 = self.d2
# could make var_pooled into attribute
var_pooled = (
(d1.sumsquares + d2.sumsquares)
/
# (d1.nobs - d1.ddof + d2.nobs - d2.ddof))
(d1.nobs - 1 + d2.nobs - 1)
)
return np.sqrt(var_pooled * (1.0 / d1.nobs + 1.0 / d2.nobs))
[docs]
def dof_satt(self):
"""degrees of freedom of Satterthwaite for unequal variance
"""
d1 = self.d1
d2 = self.d2
# this follows blindly the SPSS manual
# except I use ``_var`` which has ddof=0
sem1 = d1._var / (d1.nobs - 1)
sem2 = d2._var / (d2.nobs - 1)
semsum = sem1 + sem2
z1 = (sem1 / semsum) ** 2 / (d1.nobs - 1)
z2 = (sem2 / semsum) ** 2 / (d2.nobs - 1)
dof = 1.0 / (z1 + z2)
return dof
[docs]
def ttest_ind(self, alternative="two-sided", usevar="pooled", value=0):
"""ttest for the null hypothesis of identical means
this should also be the same as onewaygls, except for ddof differences
Parameters
----------
x1 : array_like, 1-D or 2-D
first of the two independent samples, see notes for 2-D case
x2 : array_like, 1-D or 2-D
second of the two independent samples, see notes for 2-D case
alternative : str
The alternative hypothesis, H1, has to be one of the following
'two-sided': H1: difference in means not equal to value (default)
'larger' : H1: difference in means larger than value
'smaller' : H1: difference in means smaller than value
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. If ``unequal``, then Welch ttest with Satterthwait degrees
of freedom is used
value : float
difference between the means under the Null hypothesis.
Returns
-------
tstat : float
test statistic
pvalue : float
pvalue of the t-test
df : int or float
degrees of freedom used in the t-test
Notes
-----
The result is independent of the user specified ddof.
"""
d1 = self.d1
d2 = self.d2
if usevar == "pooled":
stdm = self.std_meandiff_pooledvar
dof = d1.nobs - 1 + d2.nobs - 1
elif usevar == "unequal":
stdm = self.std_meandiff_separatevar
dof = self.dof_satt()
else:
raise ValueError('usevar can only be "pooled" or "unequal"')
tstat, pval = _tstat_generic(
d1.mean, d2.mean, stdm, dof, alternative, diff=value
)
return tstat, pval, dof
[docs]
def ztest_ind(self, alternative="two-sided", usevar="pooled", value=0):
"""z-test for the null hypothesis of identical means
Parameters
----------
x1 : array_like, 1-D or 2-D
first of the two independent samples, see notes for 2-D case
x2 : array_like, 1-D or 2-D
second of the two independent samples, see notes for 2-D case
alternative : str
The alternative hypothesis, H1, has to be one of the following
'two-sided': H1: difference in means not equal to value (default)
'larger' : H1: difference in means larger than value
'smaller' : H1: difference in means smaller than value
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. If ``unequal``, then the standard deviations of the samples may
be different.
value : float
difference between the means under the Null hypothesis.
Returns
-------
tstat : float
test statistic
pvalue : float
pvalue of the z-test
"""
d1 = self.d1
d2 = self.d2
if usevar == "pooled":
stdm = self.std_meandiff_pooledvar
elif usevar == "unequal":
stdm = self.std_meandiff_separatevar
else:
raise ValueError('usevar can only be "pooled" or "unequal"')
tstat, pval = _zstat_generic(
d1.mean, d2.mean, stdm, alternative, diff=value
)
return tstat, pval
[docs]
def tconfint_diff(
self, alpha=0.05, alternative="two-sided", usevar="pooled"
):
"""confidence interval for the difference in means
Parameters
----------
alpha : float
significance level for the confidence interval, coverage is
``1-alpha``
alternative : str
This specifies the alternative hypothesis for the test that
corresponds to the confidence interval.
The alternative hypothesis, H1, has to be one of the following :
'two-sided': H1: difference in means not equal to value (default)
'larger' : H1: difference in means larger than value
'smaller' : H1: difference in means smaller than value
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. If ``unequal``, then Welch ttest with Satterthwait degrees
of freedom is used
Returns
-------
lower, upper : floats
lower and upper limits of the confidence interval
Notes
-----
The result is independent of the user specified ddof.
"""
d1 = self.d1
d2 = self.d2
diff = d1.mean - d2.mean
if usevar == "pooled":
std_diff = self.std_meandiff_pooledvar
dof = d1.nobs - 1 + d2.nobs - 1
elif usevar == "unequal":
std_diff = self.std_meandiff_separatevar
dof = self.dof_satt()
else:
raise ValueError('usevar can only be "pooled" or "unequal"')
res = _tconfint_generic(
diff, std_diff, dof, alpha=alpha, alternative=alternative
)
return res
[docs]
def zconfint_diff(
self, alpha=0.05, alternative="two-sided", usevar="pooled"
):
"""confidence interval for the difference in means
Parameters
----------
alpha : float
significance level for the confidence interval, coverage is
``1-alpha``
alternative : str
This specifies the alternative hypothesis for the test that
corresponds to the confidence interval.
The alternative hypothesis, H1, has to be one of the following :
'two-sided': H1: difference in means not equal to value (default)
'larger' : H1: difference in means larger than value
'smaller' : H1: difference in means smaller than value
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. If ``unequal``, then Welch ttest with Satterthwait degrees
of freedom is used
Returns
-------
lower, upper : floats
lower and upper limits of the confidence interval
Notes
-----
The result is independent of the user specified ddof.
"""
d1 = self.d1
d2 = self.d2
diff = d1.mean - d2.mean
if usevar == "pooled":
std_diff = self.std_meandiff_pooledvar
elif usevar == "unequal":
std_diff = self.std_meandiff_separatevar
else:
raise ValueError('usevar can only be "pooled" or "unequal"')
res = _zconfint_generic(
diff, std_diff, alpha=alpha, alternative=alternative
)
return res
[docs]
def ttost_ind(self, low, upp, usevar="pooled"):
"""
test of equivalence for two independent samples, base on t-test
Parameters
----------
low, upp : float
equivalence interval low < m1 - m2 < upp
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. If ``unequal``, then Welch ttest with Satterthwait degrees
of freedom is used
Returns
-------
pvalue : float
pvalue of the non-equivalence test
t1, pv1 : tuple of floats
test statistic and pvalue for lower threshold test
t2, pv2 : tuple of floats
test statistic and pvalue for upper threshold test
"""
tt1 = self.ttest_ind(alternative="larger", usevar=usevar, value=low)
tt2 = self.ttest_ind(alternative="smaller", usevar=usevar, value=upp)
# TODO: remove tuple return, use same as for function tost_ind
return np.maximum(tt1[1], tt2[1]), (tt1, tt2)
[docs]
def ztost_ind(self, low, upp, usevar="pooled"):
"""
test of equivalence for two independent samples, based on z-test
Parameters
----------
low, upp : float
equivalence interval low < m1 - m2 < upp
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. If ``unequal``, then Welch ttest with Satterthwait degrees
of freedom is used
Returns
-------
pvalue : float
pvalue of the non-equivalence test
t1, pv1 : tuple of floats
test statistic and pvalue for lower threshold test
t2, pv2 : tuple of floats
test statistic and pvalue for upper threshold test
"""
tt1 = self.ztest_ind(alternative="larger", usevar=usevar, value=low)
tt2 = self.ztest_ind(alternative="smaller", usevar=usevar, value=upp)
# TODO: remove tuple return, use same as for function tost_ind
return np.maximum(tt1[1], tt2[1]), tt1, tt2
# tost.__doc__ = tost_ind.__doc__
# does not work for 2d, does not take weights into account
## def test_equal_var(self):
## """Levene test for independence
##
## """
## d1 = self.d1
## d2 = self.d2
## #rewrite this, for now just use scipy.stats
## return stats.levene(d1.data, d2.data)
[docs]
def ttest_ind(
x1,
x2,
alternative="two-sided",
usevar="pooled",
weights=(None, None),
value=0,
):
"""ttest independent sample
Convenience function that uses the classes and throws away the intermediate
results,
compared to scipy stats: drops axis option, adds alternative, usevar, and
weights option.
Parameters
----------
x1 : array_like, 1-D or 2-D
first of the two independent samples, see notes for 2-D case
x2 : array_like, 1-D or 2-D
second of the two independent samples, see notes for 2-D case
alternative : str
The alternative hypothesis, H1, has to be one of the following
* 'two-sided' (default): H1: difference in means not equal to value
* 'larger' : H1: difference in means larger than value
* 'smaller' : H1: difference in means smaller than value
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. If ``unequal``, then Welch ttest with Satterthwait degrees
of freedom is used
weights : tuple of None or ndarrays
Case weights for the two samples. For details on weights see
``DescrStatsW``
value : float
difference between the means under the Null hypothesis.
Returns
-------
tstat : float
test statistic
pvalue : float
pvalue of the t-test
df : int or float
degrees of freedom used in the t-test
"""
cm = CompareMeans(
DescrStatsW(x1, weights=weights[0], ddof=0),
DescrStatsW(x2, weights=weights[1], ddof=0),
)
tstat, pval, dof = cm.ttest_ind(
alternative=alternative, usevar=usevar, value=value
)
return tstat, pval, dof
[docs]
def ttost_ind(
x1, x2, low, upp, usevar="pooled", weights=(None, None), transform=None
):
"""test of (non-)equivalence for two independent samples
TOST: two one-sided t tests
null hypothesis: m1 - m2 < low or m1 - m2 > upp
alternative hypothesis: low < m1 - m2 < upp
where m1, m2 are the means, expected values of the two samples.
If the pvalue is smaller than a threshold, say 0.05, then we reject the
hypothesis that the difference between the two samples is larger than the
the thresholds given by low and upp.
Parameters
----------
x1 : array_like, 1-D or 2-D
first of the two independent samples, see notes for 2-D case
x2 : array_like, 1-D or 2-D
second of the two independent samples, see notes for 2-D case
low, upp : float
equivalence interval low < m1 - m2 < upp
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. If ``unequal``, then Welch ttest with Satterthwait degrees
of freedom is used
weights : tuple of None or ndarrays
Case weights for the two samples. For details on weights see
``DescrStatsW``
transform : None or function
If None (default), then the data is not transformed. Given a function,
sample data and thresholds are transformed. If transform is log, then
the equivalence interval is in ratio: low < m1 / m2 < upp
Returns
-------
pvalue : float
pvalue of the non-equivalence test
t1, pv1 : tuple of floats
test statistic and pvalue for lower threshold test
t2, pv2 : tuple of floats
test statistic and pvalue for upper threshold test
Notes
-----
The test rejects if the 2*alpha confidence interval for the difference
is contained in the ``(low, upp)`` interval.
This test works also for multi-endpoint comparisons: If d1 and d2
have the same number of columns, then each column of the data in d1 is
compared with the corresponding column in d2. This is the same as
comparing each of the corresponding columns separately. Currently no
multi-comparison correction is used. The raw p-values reported here can
be correction with the functions in ``multitest``.
"""
if transform:
if transform is np.log:
# avoid hstack in special case
x1 = transform(x1)
x2 = transform(x2)
else:
# for transforms like rankdata that will need both datasets
# concatenate works for stacking 1d and 2d arrays
xx = transform(np.concatenate((x1, x2), 0))
x1 = xx[: len(x1)]
x2 = xx[len(x1) :]
low = transform(low)
upp = transform(upp)
cm = CompareMeans(
DescrStatsW(x1, weights=weights[0], ddof=0),
DescrStatsW(x2, weights=weights[1], ddof=0),
)
pval, res = cm.ttost_ind(low, upp, usevar=usevar)
return pval, res[0], res[1]
[docs]
def ttost_paired(x1, x2, low, upp, transform=None, weights=None):
"""test of (non-)equivalence for two dependent, paired sample
TOST: two one-sided t tests
null hypothesis: md < low or md > upp
alternative hypothesis: low < md < upp
where md is the mean, expected value of the difference x1 - x2
If the pvalue is smaller than a threshold,say 0.05, then we reject the
hypothesis that the difference between the two samples is larger than the
the thresholds given by low and upp.
Parameters
----------
x1 : array_like
first of the two independent samples
x2 : array_like
second of the two independent samples
low, upp : float
equivalence interval low < mean of difference < upp
weights : None or ndarray
case weights for the two samples. For details on weights see
``DescrStatsW``
transform : None or function
If None (default), then the data is not transformed. Given a function
sample data and thresholds are transformed. If transform is log the
the equivalence interval is in ratio: low < x1 / x2 < upp
Returns
-------
pvalue : float
pvalue of the non-equivalence test
t1, pv1, df1 : tuple
test statistic, pvalue and degrees of freedom for lower threshold test
t2, pv2, df2 : tuple
test statistic, pvalue and degrees of freedom for upper threshold test
"""
if transform:
if transform is np.log:
# avoid hstack in special case
x1 = transform(x1)
x2 = transform(x2)
else:
# for transforms like rankdata that will need both datasets
# concatenate works for stacking 1d and 2d arrays
xx = transform(np.concatenate((x1, x2), 0))
x1 = xx[: len(x1)]
x2 = xx[len(x1) :]
low = transform(low)
upp = transform(upp)
dd = DescrStatsW(x1 - x2, weights=weights, ddof=0)
t1, pv1, df1 = dd.ttest_mean(low, alternative="larger")
t2, pv2, df2 = dd.ttest_mean(upp, alternative="smaller")
return np.maximum(pv1, pv2), (t1, pv1, df1), (t2, pv2, df2)
[docs]
def ztest(
x1, x2=None, value=0, alternative="two-sided", usevar="pooled", ddof=1.0
):
"""test for mean based on normal distribution, one or two samples
In the case of two samples, the samples are assumed to be independent.
Parameters
----------
x1 : array_like, 1-D or 2-D
first of the two independent samples
x2 : array_like, 1-D or 2-D
second of the two independent samples
value : float
In the one sample case, value is the mean of x1 under the Null
hypothesis.
In the two sample case, value is the difference between mean of x1 and
mean of x2 under the Null hypothesis. The test statistic is
`x1_mean - x2_mean - value`.
alternative : str
The alternative hypothesis, H1, has to be one of the following
'two-sided': H1: difference in means not equal to value (default)
'larger' : H1: difference in means larger than value
'smaller' : H1: difference in means smaller than value
usevar : str, 'pooled' or 'unequal'
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. If ``unequal``, then the standard deviation of the sample is
assumed to be different.
ddof : int
Degrees of freedom use in the calculation of the variance of the mean
estimate. In the case of comparing means this is one, however it can
be adjusted for testing other statistics (proportion, correlation)
Returns
-------
tstat : float
test statistic
pvalue : float
pvalue of the z-test
Notes
-----
usevar can be pooled or unequal in two sample case
"""
# TODO: this should delegate to CompareMeans like ttest_ind
# However that does not implement ddof
# usevar can be pooled or unequal
if usevar not in {"pooled", "unequal"}:
raise NotImplementedError('usevar can only be "pooled" or "unequal"')
x1 = np.asarray(x1)
nobs1 = x1.shape[0]
x1_mean = x1.mean(0)
x1_var = x1.var(0)
if x2 is not None:
x2 = np.asarray(x2)
nobs2 = x2.shape[0]
x2_mean = x2.mean(0)
x2_var = x2.var(0)
if usevar == "pooled":
var = nobs1 * x1_var + nobs2 * x2_var
var /= nobs1 + nobs2 - 2 * ddof
var *= 1.0 / nobs1 + 1.0 / nobs2
elif usevar == "unequal":
var = x1_var / (nobs1 - ddof) + x2_var / (nobs2 - ddof)
else:
var = x1_var / (nobs1 - ddof)
x2_mean = 0
std_diff = np.sqrt(var)
# stat = x1_mean - x2_mean - value
return _zstat_generic(x1_mean, x2_mean, std_diff, alternative, diff=value)
[docs]
def zconfint(
x1,
x2=None,
value=0,
alpha=0.05,
alternative="two-sided",
usevar="pooled",
ddof=1.0,
):
"""confidence interval based on normal distribution z-test
Parameters
----------
x1 : array_like, 1-D or 2-D
first of the two independent samples, see notes for 2-D case
x2 : array_like, 1-D or 2-D
second of the two independent samples, see notes for 2-D case
value : float
In the one sample case, value is the mean of x1 under the Null
hypothesis.
In the two sample case, value is the difference between mean of x1 and
mean of x2 under the Null hypothesis. The test statistic is
`x1_mean - x2_mean - value`.
usevar : str, 'pooled'
Currently, only 'pooled' is implemented.
If ``pooled``, then the standard deviation of the samples is assumed to be
the same. see CompareMeans.ztest_ind for different options.
ddof : int
Degrees of freedom use in the calculation of the variance of the mean
estimate. In the case of comparing means this is one, however it can
be adjusted for testing other statistics (proportion, correlation)
Notes
-----
checked only for 1 sample case
usevar not implemented, is always pooled in two sample case
``value`` shifts the confidence interval so it is centered at
`x1_mean - x2_mean - value`
See Also
--------
ztest
CompareMeans
"""
# usevar is not used, always pooled
# mostly duplicate code from ztest
if usevar != "pooled":
raise NotImplementedError('only usevar="pooled" is implemented')
x1 = np.asarray(x1)
nobs1 = x1.shape[0]
x1_mean = x1.mean(0)
x1_var = x1.var(0)
if x2 is not None:
x2 = np.asarray(x2)
nobs2 = x2.shape[0]
x2_mean = x2.mean(0)
x2_var = x2.var(0)
var_pooled = nobs1 * x1_var + nobs2 * x2_var
var_pooled /= nobs1 + nobs2 - 2 * ddof
var_pooled *= 1.0 / nobs1 + 1.0 / nobs2
else:
var_pooled = x1_var / (nobs1 - ddof)
x2_mean = 0
std_diff = np.sqrt(var_pooled)
ci = _zconfint_generic(
x1_mean - x2_mean - value, std_diff, alpha, alternative
)
return ci
[docs]
def ztost(x1, low, upp, x2=None, usevar="pooled", ddof=1.0):
"""Equivalence test based on normal distribution
Parameters
----------
x1 : array_like
one sample or first sample for 2 independent samples
low, upp : float
equivalence interval low < m1 - m2 < upp
x1 : array_like or None
second sample for 2 independent samples test. If None, then a
one-sample test is performed.
usevar : str, 'pooled'
If `pooled`, then the standard deviation of the samples is assumed to be
the same. Only `pooled` is currently implemented.
Returns
-------
pvalue : float
pvalue of the non-equivalence test
t1, pv1 : tuple of floats
test statistic and pvalue for lower threshold test
t2, pv2 : tuple of floats
test statistic and pvalue for upper threshold test
Notes
-----
checked only for 1 sample case
"""
tt1 = ztest(
x1, x2, alternative="larger", usevar=usevar, value=low, ddof=ddof
)
tt2 = ztest(
x1, x2, alternative="smaller", usevar=usevar, value=upp, ddof=ddof
)
return (
np.maximum(tt1[1], tt2[1]),
tt1,
tt2,
)
Last update:
Nov 14, 2024