Coverage for /usr/lib/python3/dist-packages/scipy/stats/_hypotests.py: 12%
477 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from collections import namedtuple
2from dataclasses import dataclass
3from math import comb
4import numpy as np
5import warnings
6from itertools import combinations
7import scipy.stats
8from scipy.optimize import shgo
9from . import distributions
10from ._common import ConfidenceInterval
11from ._continuous_distns import chi2, norm
12from scipy.special import gamma, kv, gammaln
13from scipy.fft import ifft
14from ._stats_pythran import _a_ij_Aij_Dij2
15from ._stats_pythran import (
16 _concordant_pairs as _P, _discordant_pairs as _Q
17)
18from scipy.stats import _stats_py
20__all__ = ['epps_singleton_2samp', 'cramervonmises', 'somersd',
21 'barnard_exact', 'boschloo_exact', 'cramervonmises_2samp',
22 'tukey_hsd', 'poisson_means_test']
24Epps_Singleton_2sampResult = namedtuple('Epps_Singleton_2sampResult',
25 ('statistic', 'pvalue'))
28def epps_singleton_2samp(x, y, t=(0.4, 0.8)):
29 """Compute the Epps-Singleton (ES) test statistic.
31 Test the null hypothesis that two samples have the same underlying
32 probability distribution.
34 Parameters
35 ----------
36 x, y : array-like
37 The two samples of observations to be tested. Input must not have more
38 than one dimension. Samples can have different lengths.
39 t : array-like, optional
40 The points (t1, ..., tn) where the empirical characteristic function is
41 to be evaluated. It should be positive distinct numbers. The default
42 value (0.4, 0.8) is proposed in [1]_. Input must not have more than
43 one dimension.
45 Returns
46 -------
47 statistic : float
48 The test statistic.
49 pvalue : float
50 The associated p-value based on the asymptotic chi2-distribution.
52 See Also
53 --------
54 ks_2samp, anderson_ksamp
56 Notes
57 -----
58 Testing whether two samples are generated by the same underlying
59 distribution is a classical question in statistics. A widely used test is
60 the Kolmogorov-Smirnov (KS) test which relies on the empirical
61 distribution function. Epps and Singleton introduce a test based on the
62 empirical characteristic function in [1]_.
64 One advantage of the ES test compared to the KS test is that is does
65 not assume a continuous distribution. In [1]_, the authors conclude
66 that the test also has a higher power than the KS test in many
67 examples. They recommend the use of the ES test for discrete samples as
68 well as continuous samples with at least 25 observations each, whereas
69 `anderson_ksamp` is recommended for smaller sample sizes in the
70 continuous case.
72 The p-value is computed from the asymptotic distribution of the test
73 statistic which follows a `chi2` distribution. If the sample size of both
74 `x` and `y` is below 25, the small sample correction proposed in [1]_ is
75 applied to the test statistic.
77 The default values of `t` are determined in [1]_ by considering
78 various distributions and finding good values that lead to a high power
79 of the test in general. Table III in [1]_ gives the optimal values for
80 the distributions tested in that study. The values of `t` are scaled by
81 the semi-interquartile range in the implementation, see [1]_.
83 References
84 ----------
85 .. [1] T. W. Epps and K. J. Singleton, "An omnibus test for the two-sample
86 problem using the empirical characteristic function", Journal of
87 Statistical Computation and Simulation 26, p. 177--203, 1986.
89 .. [2] S. J. Goerg and J. Kaiser, "Nonparametric testing of distributions
90 - the Epps-Singleton two-sample test using the empirical characteristic
91 function", The Stata Journal 9(3), p. 454--465, 2009.
93 """
94 x, y, t = np.asarray(x), np.asarray(y), np.asarray(t)
95 # check if x and y are valid inputs
96 if x.ndim > 1:
97 raise ValueError(f'x must be 1d, but x.ndim equals {x.ndim}.')
98 if y.ndim > 1:
99 raise ValueError(f'y must be 1d, but y.ndim equals {y.ndim}.')
100 nx, ny = len(x), len(y)
101 if (nx < 5) or (ny < 5):
102 raise ValueError('x and y should have at least 5 elements, but len(x) '
103 '= {} and len(y) = {}.'.format(nx, ny))
104 if not np.isfinite(x).all():
105 raise ValueError('x must not contain nonfinite values.')
106 if not np.isfinite(y).all():
107 raise ValueError('y must not contain nonfinite values.')
108 n = nx + ny
110 # check if t is valid
111 if t.ndim > 1:
112 raise ValueError(f't must be 1d, but t.ndim equals {t.ndim}.')
113 if np.less_equal(t, 0).any():
114 raise ValueError('t must contain positive elements only.')
116 # rescale t with semi-iqr as proposed in [1]; import iqr here to avoid
117 # circular import
118 from scipy.stats import iqr
119 sigma = iqr(np.hstack((x, y))) / 2
120 ts = np.reshape(t, (-1, 1)) / sigma
122 # covariance estimation of ES test
123 gx = np.vstack((np.cos(ts*x), np.sin(ts*x))).T # shape = (nx, 2*len(t))
124 gy = np.vstack((np.cos(ts*y), np.sin(ts*y))).T
125 cov_x = np.cov(gx.T, bias=True) # the test uses biased cov-estimate
126 cov_y = np.cov(gy.T, bias=True)
127 est_cov = (n/nx)*cov_x + (n/ny)*cov_y
128 est_cov_inv = np.linalg.pinv(est_cov)
129 r = np.linalg.matrix_rank(est_cov_inv)
130 if r < 2*len(t):
131 warnings.warn('Estimated covariance matrix does not have full rank. '
132 'This indicates a bad choice of the input t and the '
133 'test might not be consistent.') # see p. 183 in [1]_
135 # compute test statistic w distributed asympt. as chisquare with df=r
136 g_diff = np.mean(gx, axis=0) - np.mean(gy, axis=0)
137 w = n*np.dot(g_diff.T, np.dot(est_cov_inv, g_diff))
139 # apply small-sample correction
140 if (max(nx, ny) < 25):
141 corr = 1.0/(1.0 + n**(-0.45) + 10.1*(nx**(-1.7) + ny**(-1.7)))
142 w = corr * w
144 p = chi2.sf(w, r)
146 return Epps_Singleton_2sampResult(w, p)
149def poisson_means_test(k1, n1, k2, n2, *, diff=0, alternative='two-sided'):
150 r"""
151 Performs the Poisson means test, AKA the "E-test".
153 This is a test of the null hypothesis that the difference between means of
154 two Poisson distributions is `diff`. The samples are provided as the
155 number of events `k1` and `k2` observed within measurement intervals
156 (e.g. of time, space, number of observations) of sizes `n1` and `n2`.
158 Parameters
159 ----------
160 k1 : int
161 Number of events observed from distribution 1.
162 n1: float
163 Size of sample from distribution 1.
164 k2 : int
165 Number of events observed from distribution 2.
166 n2 : float
167 Size of sample from distribution 2.
168 diff : float, default=0
169 The hypothesized difference in means between the distributions
170 underlying the samples.
171 alternative : {'two-sided', 'less', 'greater'}, optional
172 Defines the alternative hypothesis.
173 The following options are available (default is 'two-sided'):
175 * 'two-sided': the difference between distribution means is not
176 equal to `diff`
177 * 'less': the difference between distribution means is less than
178 `diff`
179 * 'greater': the difference between distribution means is greater
180 than `diff`
182 Returns
183 -------
184 statistic : float
185 The test statistic (see [1]_ equation 3.3).
186 pvalue : float
187 The probability of achieving such an extreme value of the test
188 statistic under the null hypothesis.
190 Notes
191 -----
193 Let:
195 .. math:: X_1 \sim \mbox{Poisson}(\mathtt{n1}\lambda_1)
197 be a random variable independent of
199 .. math:: X_2 \sim \mbox{Poisson}(\mathtt{n2}\lambda_2)
201 and let ``k1`` and ``k2`` be the observed values of :math:`X_1`
202 and :math:`X_2`, respectively. Then `poisson_means_test` uses the number
203 of observed events ``k1`` and ``k2`` from samples of size ``n1`` and
204 ``n2``, respectively, to test the null hypothesis that
206 .. math::
207 H_0: \lambda_1 - \lambda_2 = \mathtt{diff}
209 A benefit of the E-test is that it has good power for small sample sizes,
210 which can reduce sampling costs [1]_. It has been evaluated and determined
211 to be more powerful than the comparable C-test, sometimes referred to as
212 the Poisson exact test.
214 References
215 ----------
216 .. [1] Krishnamoorthy, K., & Thomson, J. (2004). A more powerful test for
217 comparing two Poisson means. Journal of Statistical Planning and
218 Inference, 119(1), 23-35.
220 .. [2] Przyborowski, J., & Wilenski, H. (1940). Homogeneity of results in
221 testing samples from Poisson series: With an application to testing
222 clover seed for dodder. Biometrika, 31(3/4), 313-323.
224 Examples
225 --------
227 Suppose that a gardener wishes to test the number of dodder (weed) seeds
228 in a sack of clover seeds that they buy from a seed company. It has
229 previously been established that the number of dodder seeds in clover
230 follows the Poisson distribution.
232 A 100 gram sample is drawn from the sack before being shipped to the
233 gardener. The sample is analyzed, and it is found to contain no dodder
234 seeds; that is, `k1` is 0. However, upon arrival, the gardener draws
235 another 100 gram sample from the sack. This time, three dodder seeds are
236 found in the sample; that is, `k2` is 3. The gardener would like to
237 know if the difference is significant and not due to chance. The
238 null hypothesis is that the difference between the two samples is merely
239 due to chance, or that :math:`\lambda_1 - \lambda_2 = \mathtt{diff}`
240 where :math:`\mathtt{diff} = 0`. The alternative hypothesis is that the
241 difference is not due to chance, or :math:`\lambda_1 - \lambda_2 \ne 0`.
242 The gardener selects a significance level of 5% to reject the null
243 hypothesis in favor of the alternative [2]_.
245 >>> import scipy.stats as stats
246 >>> res = stats.poisson_means_test(0, 100, 3, 100)
247 >>> res.statistic, res.pvalue
248 (-1.7320508075688772, 0.08837900929018157)
250 The p-value is .088, indicating a near 9% chance of observing a value of
251 the test statistic under the null hypothesis. This exceeds 5%, so the
252 gardener does not reject the null hypothesis as the difference cannot be
253 regarded as significant at this level.
254 """
256 _poisson_means_test_iv(k1, n1, k2, n2, diff, alternative)
258 # "for a given k_1 and k_2, an estimate of \lambda_2 is given by" [1] (3.4)
259 lmbd_hat2 = ((k1 + k2) / (n1 + n2) - diff * n1 / (n1 + n2))
261 # "\hat{\lambda_{2k}} may be less than or equal to zero ... and in this
262 # case the null hypothesis cannot be rejected ... [and] it is not necessary
263 # to compute the p-value". [1] page 26 below eq. (3.6).
264 if lmbd_hat2 <= 0:
265 return _stats_py.SignificanceResult(0, 1)
267 # The unbiased variance estimate [1] (3.2)
268 var = k1 / (n1 ** 2) + k2 / (n2 ** 2)
270 # The _observed_ pivot statistic from the input. It follows the
271 # unnumbered equation following equation (3.3) This is used later in
272 # comparison with the computed pivot statistics in an indicator function.
273 t_k1k2 = (k1 / n1 - k2 / n2 - diff) / np.sqrt(var)
275 # Equation (3.5) of [1] is lengthy, so it is broken into several parts,
276 # beginning here. Note that the probability mass function of poisson is
277 # exp^(-\mu)*\mu^k/k!, so and this is called with shape \mu, here noted
278 # here as nlmbd_hat*. The strategy for evaluating the double summation in
279 # (3.5) is to create two arrays of the values of the two products inside
280 # the summation and then broadcast them together into a matrix, and then
281 # sum across the entire matrix.
283 # Compute constants (as seen in the first and second separated products in
284 # (3.5).). (This is the shape (\mu) parameter of the poisson distribution.)
285 nlmbd_hat1 = n1 * (lmbd_hat2 + diff)
286 nlmbd_hat2 = n2 * lmbd_hat2
288 # Determine summation bounds for tail ends of distribution rather than
289 # summing to infinity. `x1*` is for the outer sum and `x2*` is the inner
290 # sum.
291 x1_lb, x1_ub = distributions.poisson.ppf([1e-10, 1 - 1e-16], nlmbd_hat1)
292 x2_lb, x2_ub = distributions.poisson.ppf([1e-10, 1 - 1e-16], nlmbd_hat2)
294 # Construct arrays to function as the x_1 and x_2 counters on the summation
295 # in (3.5). `x1` is in columns and `x2` is in rows to allow for
296 # broadcasting.
297 x1 = np.arange(x1_lb, x1_ub + 1)
298 x2 = np.arange(x2_lb, x2_ub + 1)[:, None]
300 # These are the two products in equation (3.5) with `prob_x1` being the
301 # first (left side) and `prob_x2` being the second (right side). (To
302 # make as clear as possible: the 1st contains a "+ d" term, the 2nd does
303 # not.)
304 prob_x1 = distributions.poisson.pmf(x1, nlmbd_hat1)
305 prob_x2 = distributions.poisson.pmf(x2, nlmbd_hat2)
307 # compute constants for use in the "pivot statistic" per the
308 # unnumbered equation following (3.3).
309 lmbd_x1 = x1 / n1
310 lmbd_x2 = x2 / n2
311 lmbds_diff = lmbd_x1 - lmbd_x2 - diff
312 var_x1x2 = lmbd_x1 / n1 + lmbd_x2 / n2
314 # This is the 'pivot statistic' for use in the indicator of the summation
315 # (left side of "I[.]").
316 with np.errstate(invalid='ignore', divide='ignore'):
317 t_x1x2 = lmbds_diff / np.sqrt(var_x1x2)
319 # `[indicator]` implements the "I[.] ... the indicator function" per
320 # the paragraph following equation (3.5).
321 if alternative == 'two-sided':
322 indicator = np.abs(t_x1x2) >= np.abs(t_k1k2)
323 elif alternative == 'less':
324 indicator = t_x1x2 <= t_k1k2
325 else:
326 indicator = t_x1x2 >= t_k1k2
328 # Multiply all combinations of the products together, exclude terms
329 # based on the `indicator` and then sum. (3.5)
330 pvalue = np.sum((prob_x1 * prob_x2)[indicator])
331 return _stats_py.SignificanceResult(t_k1k2, pvalue)
334def _poisson_means_test_iv(k1, n1, k2, n2, diff, alternative):
335 # """check for valid types and values of input to `poisson_mean_test`."""
336 if k1 != int(k1) or k2 != int(k2):
337 raise TypeError('`k1` and `k2` must be integers.')
339 count_err = '`k1` and `k2` must be greater than or equal to 0.'
340 if k1 < 0 or k2 < 0:
341 raise ValueError(count_err)
343 if n1 <= 0 or n2 <= 0:
344 raise ValueError('`n1` and `n2` must be greater than 0.')
346 if diff < 0:
347 raise ValueError('diff must be greater than or equal to 0.')
349 alternatives = {'two-sided', 'less', 'greater'}
350 if alternative.lower() not in alternatives:
351 raise ValueError(f"Alternative must be one of '{alternatives}'.")
354class CramerVonMisesResult:
355 def __init__(self, statistic, pvalue):
356 self.statistic = statistic
357 self.pvalue = pvalue
359 def __repr__(self):
360 return (f"{self.__class__.__name__}(statistic={self.statistic}, "
361 f"pvalue={self.pvalue})")
364def _psi1_mod(x):
365 """
366 psi1 is defined in equation 1.10 in Csörgő, S. and Faraway, J. (1996).
367 This implements a modified version by excluding the term V(x) / 12
368 (here: _cdf_cvm_inf(x) / 12) to avoid evaluating _cdf_cvm_inf(x)
369 twice in _cdf_cvm.
371 Implementation based on MAPLE code of Julian Faraway and R code of the
372 function pCvM in the package goftest (v1.1.1), permission granted
373 by Adrian Baddeley. Main difference in the implementation: the code
374 here keeps adding terms of the series until the terms are small enough.
375 """
377 def _ed2(y):
378 z = y**2 / 4
379 b = kv(1/4, z) + kv(3/4, z)
380 return np.exp(-z) * (y/2)**(3/2) * b / np.sqrt(np.pi)
382 def _ed3(y):
383 z = y**2 / 4
384 c = np.exp(-z) / np.sqrt(np.pi)
385 return c * (y/2)**(5/2) * (2*kv(1/4, z) + 3*kv(3/4, z) - kv(5/4, z))
387 def _Ak(k, x):
388 m = 2*k + 1
389 sx = 2 * np.sqrt(x)
390 y1 = x**(3/4)
391 y2 = x**(5/4)
393 e1 = m * gamma(k + 1/2) * _ed2((4 * k + 3)/sx) / (9 * y1)
394 e2 = gamma(k + 1/2) * _ed3((4 * k + 1) / sx) / (72 * y2)
395 e3 = 2 * (m + 2) * gamma(k + 3/2) * _ed3((4 * k + 5) / sx) / (12 * y2)
396 e4 = 7 * m * gamma(k + 1/2) * _ed2((4 * k + 1) / sx) / (144 * y1)
397 e5 = 7 * m * gamma(k + 1/2) * _ed2((4 * k + 5) / sx) / (144 * y1)
399 return e1 + e2 + e3 + e4 + e5
401 x = np.asarray(x)
402 tot = np.zeros_like(x, dtype='float')
403 cond = np.ones_like(x, dtype='bool')
404 k = 0
405 while np.any(cond):
406 z = -_Ak(k, x[cond]) / (np.pi * gamma(k + 1))
407 tot[cond] = tot[cond] + z
408 cond[cond] = np.abs(z) >= 1e-7
409 k += 1
411 return tot
414def _cdf_cvm_inf(x):
415 """
416 Calculate the cdf of the Cramér-von Mises statistic (infinite sample size).
418 See equation 1.2 in Csörgő, S. and Faraway, J. (1996).
420 Implementation based on MAPLE code of Julian Faraway and R code of the
421 function pCvM in the package goftest (v1.1.1), permission granted
422 by Adrian Baddeley. Main difference in the implementation: the code
423 here keeps adding terms of the series until the terms are small enough.
425 The function is not expected to be accurate for large values of x, say
426 x > 4, when the cdf is very close to 1.
427 """
428 x = np.asarray(x)
430 def term(x, k):
431 # this expression can be found in [2], second line of (1.3)
432 u = np.exp(gammaln(k + 0.5) - gammaln(k+1)) / (np.pi**1.5 * np.sqrt(x))
433 y = 4*k + 1
434 q = y**2 / (16*x)
435 b = kv(0.25, q)
436 return u * np.sqrt(y) * np.exp(-q) * b
438 tot = np.zeros_like(x, dtype='float')
439 cond = np.ones_like(x, dtype='bool')
440 k = 0
441 while np.any(cond):
442 z = term(x[cond], k)
443 tot[cond] = tot[cond] + z
444 cond[cond] = np.abs(z) >= 1e-7
445 k += 1
447 return tot
450def _cdf_cvm(x, n=None):
451 """
452 Calculate the cdf of the Cramér-von Mises statistic for a finite sample
453 size n. If N is None, use the asymptotic cdf (n=inf).
455 See equation 1.8 in Csörgő, S. and Faraway, J. (1996) for finite samples,
456 1.2 for the asymptotic cdf.
458 The function is not expected to be accurate for large values of x, say
459 x > 2, when the cdf is very close to 1 and it might return values > 1
460 in that case, e.g. _cdf_cvm(2.0, 12) = 1.0000027556716846. Moreover, it
461 is not accurate for small values of n, especially close to the bounds of
462 the distribution's domain, [1/(12*n), n/3], where the value jumps to 0
463 and 1, respectively. These are limitations of the approximation by Csörgő
464 and Faraway (1996) implemented in this function.
465 """
466 x = np.asarray(x)
467 if n is None:
468 y = _cdf_cvm_inf(x)
469 else:
470 # support of the test statistic is [12/n, n/3], see 1.1 in [2]
471 y = np.zeros_like(x, dtype='float')
472 sup = (1./(12*n) < x) & (x < n/3.)
473 # note: _psi1_mod does not include the term _cdf_cvm_inf(x) / 12
474 # therefore, we need to add it here
475 y[sup] = _cdf_cvm_inf(x[sup]) * (1 + 1./(12*n)) + _psi1_mod(x[sup]) / n
476 y[x >= n/3] = 1
478 if y.ndim == 0:
479 return y[()]
480 return y
483def cramervonmises(rvs, cdf, args=()):
484 """Perform the one-sample Cramér-von Mises test for goodness of fit.
486 This performs a test of the goodness of fit of a cumulative distribution
487 function (cdf) :math:`F` compared to the empirical distribution function
488 :math:`F_n` of observed random variates :math:`X_1, ..., X_n` that are
489 assumed to be independent and identically distributed ([1]_).
490 The null hypothesis is that the :math:`X_i` have cumulative distribution
491 :math:`F`.
493 Parameters
494 ----------
495 rvs : array_like
496 A 1-D array of observed values of the random variables :math:`X_i`.
497 cdf : str or callable
498 The cumulative distribution function :math:`F` to test the
499 observations against. If a string, it should be the name of a
500 distribution in `scipy.stats`. If a callable, that callable is used
501 to calculate the cdf: ``cdf(x, *args) -> float``.
502 args : tuple, optional
503 Distribution parameters. These are assumed to be known; see Notes.
505 Returns
506 -------
507 res : object with attributes
508 statistic : float
509 Cramér-von Mises statistic.
510 pvalue : float
511 The p-value.
513 See Also
514 --------
515 kstest, cramervonmises_2samp
517 Notes
518 -----
519 .. versionadded:: 1.6.0
521 The p-value relies on the approximation given by equation 1.8 in [2]_.
522 It is important to keep in mind that the p-value is only accurate if
523 one tests a simple hypothesis, i.e. the parameters of the reference
524 distribution are known. If the parameters are estimated from the data
525 (composite hypothesis), the computed p-value is not reliable.
527 References
528 ----------
529 .. [1] Cramér-von Mises criterion, Wikipedia,
530 https://en.wikipedia.org/wiki/Cram%C3%A9r%E2%80%93von_Mises_criterion
531 .. [2] Csörgő, S. and Faraway, J. (1996). The Exact and Asymptotic
532 Distribution of Cramér-von Mises Statistics. Journal of the
533 Royal Statistical Society, pp. 221-234.
535 Examples
536 --------
538 Suppose we wish to test whether data generated by ``scipy.stats.norm.rvs``
539 were, in fact, drawn from the standard normal distribution. We choose a
540 significance level of ``alpha=0.05``.
542 >>> import numpy as np
543 >>> from scipy import stats
544 >>> rng = np.random.default_rng(165417232101553420507139617764912913465)
545 >>> x = stats.norm.rvs(size=500, random_state=rng)
546 >>> res = stats.cramervonmises(x, 'norm')
547 >>> res.statistic, res.pvalue
548 (0.1072085112565724, 0.5508482238203407)
550 The p-value exceeds our chosen significance level, so we do not
551 reject the null hypothesis that the observed sample is drawn from the
552 standard normal distribution.
554 Now suppose we wish to check whether the same samples shifted by 2.1 is
555 consistent with being drawn from a normal distribution with a mean of 2.
557 >>> y = x + 2.1
558 >>> res = stats.cramervonmises(y, 'norm', args=(2,))
559 >>> res.statistic, res.pvalue
560 (0.8364446265294695, 0.00596286797008283)
562 Here we have used the `args` keyword to specify the mean (``loc``)
563 of the normal distribution to test the data against. This is equivalent
564 to the following, in which we create a frozen normal distribution with
565 mean 2.1, then pass its ``cdf`` method as an argument.
567 >>> frozen_dist = stats.norm(loc=2)
568 >>> res = stats.cramervonmises(y, frozen_dist.cdf)
569 >>> res.statistic, res.pvalue
570 (0.8364446265294695, 0.00596286797008283)
572 In either case, we would reject the null hypothesis that the observed
573 sample is drawn from a normal distribution with a mean of 2 (and default
574 variance of 1) because the p-value is less than our chosen
575 significance level.
577 """
578 if isinstance(cdf, str):
579 cdf = getattr(distributions, cdf).cdf
581 vals = np.sort(np.asarray(rvs))
583 if vals.size <= 1:
584 raise ValueError('The sample must contain at least two observations.')
585 if vals.ndim > 1:
586 raise ValueError('The sample must be one-dimensional.')
588 n = len(vals)
589 cdfvals = cdf(vals, *args)
591 u = (2*np.arange(1, n+1) - 1)/(2*n)
592 w = 1/(12*n) + np.sum((u - cdfvals)**2)
594 # avoid small negative values that can occur due to the approximation
595 p = max(0, 1. - _cdf_cvm(w, n))
597 return CramerVonMisesResult(statistic=w, pvalue=p)
600def _get_wilcoxon_distr(n):
601 """
602 Distribution of probability of the Wilcoxon ranksum statistic r_plus (sum
603 of ranks of positive differences).
604 Returns an array with the probabilities of all the possible ranks
605 r = 0, ..., n*(n+1)/2
606 """
607 c = np.ones(1, dtype=np.double)
608 for k in range(1, n + 1):
609 prev_c = c
610 c = np.zeros(k * (k + 1) // 2 + 1, dtype=np.double)
611 m = len(prev_c)
612 c[:m] = prev_c * 0.5
613 c[-m:] += prev_c * 0.5
614 return c
617def _get_wilcoxon_distr2(n):
618 """
619 Distribution of probability of the Wilcoxon ranksum statistic r_plus (sum
620 of ranks of positive differences).
621 Returns an array with the probabilities of all the possible ranks
622 r = 0, ..., n*(n+1)/2
623 This is a slower reference function
624 References
625 ----------
626 .. [1] 1. Harris T, Hardin JW. Exact Wilcoxon Signed-Rank and Wilcoxon
627 Mann-Whitney Ranksum Tests. The Stata Journal. 2013;13(2):337-343.
628 """
629 ai = np.arange(1, n+1)[:, None]
630 t = n*(n+1)/2
631 q = 2*t
632 j = np.arange(q)
633 theta = 2*np.pi/q*j
634 phi_sp = np.prod(np.cos(theta*ai), axis=0)
635 phi_s = np.exp(1j*theta*t) * phi_sp
636 p = np.real(ifft(phi_s))
637 res = np.zeros(int(t)+1)
638 res[:-1:] = p[::2]
639 res[0] /= 2
640 res[-1] = res[0]
641 return res
644def _tau_b(A):
645 """Calculate Kendall's tau-b and p-value from contingency table."""
646 # See [2] 2.2 and 4.2
648 # contingency table must be truly 2D
649 if A.shape[0] == 1 or A.shape[1] == 1:
650 return np.nan, np.nan
652 NA = A.sum()
653 PA = _P(A)
654 QA = _Q(A)
655 Sri2 = (A.sum(axis=1)**2).sum()
656 Scj2 = (A.sum(axis=0)**2).sum()
657 denominator = (NA**2 - Sri2)*(NA**2 - Scj2)
659 tau = (PA-QA)/(denominator)**0.5
661 numerator = 4*(_a_ij_Aij_Dij2(A) - (PA - QA)**2 / NA)
662 s02_tau_b = numerator/denominator
663 if s02_tau_b == 0: # Avoid divide by zero
664 return tau, 0
665 Z = tau/s02_tau_b**0.5
666 p = 2*norm.sf(abs(Z)) # 2-sided p-value
668 return tau, p
671def _somers_d(A, alternative='two-sided'):
672 """Calculate Somers' D and p-value from contingency table."""
673 # See [3] page 1740
675 # contingency table must be truly 2D
676 if A.shape[0] <= 1 or A.shape[1] <= 1:
677 return np.nan, np.nan
679 NA = A.sum()
680 NA2 = NA**2
681 PA = _P(A)
682 QA = _Q(A)
683 Sri2 = (A.sum(axis=1)**2).sum()
685 d = (PA - QA)/(NA2 - Sri2)
687 S = _a_ij_Aij_Dij2(A) - (PA-QA)**2/NA
689 with np.errstate(divide='ignore'):
690 Z = (PA - QA)/(4*(S))**0.5
692 _, p = scipy.stats._stats_py._normtest_finish(Z, alternative)
694 return d, p
697@dataclass
698class SomersDResult:
699 statistic: float
700 pvalue: float
701 table: np.ndarray
704def somersd(x, y=None, alternative='two-sided'):
705 r"""Calculates Somers' D, an asymmetric measure of ordinal association.
707 Like Kendall's :math:`\tau`, Somers' :math:`D` is a measure of the
708 correspondence between two rankings. Both statistics consider the
709 difference between the number of concordant and discordant pairs in two
710 rankings :math:`X` and :math:`Y`, and both are normalized such that values
711 close to 1 indicate strong agreement and values close to -1 indicate
712 strong disagreement. They differ in how they are normalized. To show the
713 relationship, Somers' :math:`D` can be defined in terms of Kendall's
714 :math:`\tau_a`:
716 .. math::
717 D(Y|X) = \frac{\tau_a(X, Y)}{\tau_a(X, X)}
719 Suppose the first ranking :math:`X` has :math:`r` distinct ranks and the
720 second ranking :math:`Y` has :math:`s` distinct ranks. These two lists of
721 :math:`n` rankings can also be viewed as an :math:`r \times s` contingency
722 table in which element :math:`i, j` is the number of rank pairs with rank
723 :math:`i` in ranking :math:`X` and rank :math:`j` in ranking :math:`Y`.
724 Accordingly, `somersd` also allows the input data to be supplied as a
725 single, 2D contingency table instead of as two separate, 1D rankings.
727 Note that the definition of Somers' :math:`D` is asymmetric: in general,
728 :math:`D(Y|X) \neq D(X|Y)`. ``somersd(x, y)`` calculates Somers'
729 :math:`D(Y|X)`: the "row" variable :math:`X` is treated as an independent
730 variable, and the "column" variable :math:`Y` is dependent. For Somers'
731 :math:`D(X|Y)`, swap the input lists or transpose the input table.
733 Parameters
734 ----------
735 x : array_like
736 1D array of rankings, treated as the (row) independent variable.
737 Alternatively, a 2D contingency table.
738 y : array_like, optional
739 If `x` is a 1D array of rankings, `y` is a 1D array of rankings of the
740 same length, treated as the (column) dependent variable.
741 If `x` is 2D, `y` is ignored.
742 alternative : {'two-sided', 'less', 'greater'}, optional
743 Defines the alternative hypothesis. Default is 'two-sided'.
744 The following options are available:
745 * 'two-sided': the rank correlation is nonzero
746 * 'less': the rank correlation is negative (less than zero)
747 * 'greater': the rank correlation is positive (greater than zero)
749 Returns
750 -------
751 res : SomersDResult
752 A `SomersDResult` object with the following fields:
754 statistic : float
755 The Somers' :math:`D` statistic.
756 pvalue : float
757 The p-value for a hypothesis test whose null
758 hypothesis is an absence of association, :math:`D=0`.
759 See notes for more information.
760 table : 2D array
761 The contingency table formed from rankings `x` and `y` (or the
762 provided contingency table, if `x` is a 2D array)
764 See Also
765 --------
766 kendalltau : Calculates Kendall's tau, another correlation measure.
767 weightedtau : Computes a weighted version of Kendall's tau.
768 spearmanr : Calculates a Spearman rank-order correlation coefficient.
769 pearsonr : Calculates a Pearson correlation coefficient.
771 Notes
772 -----
773 This function follows the contingency table approach of [2]_ and
774 [3]_. *p*-values are computed based on an asymptotic approximation of
775 the test statistic distribution under the null hypothesis :math:`D=0`.
777 Theoretically, hypothesis tests based on Kendall's :math:`tau` and Somers'
778 :math:`D` should be identical.
779 However, the *p*-values returned by `kendalltau` are based
780 on the null hypothesis of *independence* between :math:`X` and :math:`Y`
781 (i.e. the population from which pairs in :math:`X` and :math:`Y` are
782 sampled contains equal numbers of all possible pairs), which is more
783 specific than the null hypothesis :math:`D=0` used here. If the null
784 hypothesis of independence is desired, it is acceptable to use the
785 *p*-value returned by `kendalltau` with the statistic returned by
786 `somersd` and vice versa. For more information, see [2]_.
788 Contingency tables are formatted according to the convention used by
789 SAS and R: the first ranking supplied (``x``) is the "row" variable, and
790 the second ranking supplied (``y``) is the "column" variable. This is
791 opposite the convention of Somers' original paper [1]_.
793 References
794 ----------
795 .. [1] Robert H. Somers, "A New Asymmetric Measure of Association for
796 Ordinal Variables", *American Sociological Review*, Vol. 27, No. 6,
797 pp. 799--811, 1962.
799 .. [2] Morton B. Brown and Jacqueline K. Benedetti, "Sampling Behavior of
800 Tests for Correlation in Two-Way Contingency Tables", *Journal of
801 the American Statistical Association* Vol. 72, No. 358, pp.
802 309--315, 1977.
804 .. [3] SAS Institute, Inc., "The FREQ Procedure (Book Excerpt)",
805 *SAS/STAT 9.2 User's Guide, Second Edition*, SAS Publishing, 2009.
807 .. [4] Laerd Statistics, "Somers' d using SPSS Statistics", *SPSS
808 Statistics Tutorials and Statistical Guides*,
809 https://statistics.laerd.com/spss-tutorials/somers-d-using-spss-statistics.php,
810 Accessed July 31, 2020.
812 Examples
813 --------
814 We calculate Somers' D for the example given in [4]_, in which a hotel
815 chain owner seeks to determine the association between hotel room
816 cleanliness and customer satisfaction. The independent variable, hotel
817 room cleanliness, is ranked on an ordinal scale: "below average (1)",
818 "average (2)", or "above average (3)". The dependent variable, customer
819 satisfaction, is ranked on a second scale: "very dissatisfied (1)",
820 "moderately dissatisfied (2)", "neither dissatisfied nor satisfied (3)",
821 "moderately satisfied (4)", or "very satisfied (5)". 189 customers
822 respond to the survey, and the results are cast into a contingency table
823 with the hotel room cleanliness as the "row" variable and customer
824 satisfaction as the "column" variable.
826 +-----+-----+-----+-----+-----+-----+
827 | | (1) | (2) | (3) | (4) | (5) |
828 +=====+=====+=====+=====+=====+=====+
829 | (1) | 27 | 25 | 14 | 7 | 0 |
830 +-----+-----+-----+-----+-----+-----+
831 | (2) | 7 | 14 | 18 | 35 | 12 |
832 +-----+-----+-----+-----+-----+-----+
833 | (3) | 1 | 3 | 2 | 7 | 17 |
834 +-----+-----+-----+-----+-----+-----+
836 For example, 27 customers assigned their room a cleanliness ranking of
837 "below average (1)" and a corresponding satisfaction of "very
838 dissatisfied (1)". We perform the analysis as follows.
840 >>> from scipy.stats import somersd
841 >>> table = [[27, 25, 14, 7, 0], [7, 14, 18, 35, 12], [1, 3, 2, 7, 17]]
842 >>> res = somersd(table)
843 >>> res.statistic
844 0.6032766111513396
845 >>> res.pvalue
846 1.0007091191074533e-27
848 The value of the Somers' D statistic is approximately 0.6, indicating
849 a positive correlation between room cleanliness and customer satisfaction
850 in the sample.
851 The *p*-value is very small, indicating a very small probability of
852 observing such an extreme value of the statistic under the null
853 hypothesis that the statistic of the entire population (from which
854 our sample of 189 customers is drawn) is zero. This supports the
855 alternative hypothesis that the true value of Somers' D for the population
856 is nonzero.
858 """
859 x, y = np.array(x), np.array(y)
860 if x.ndim == 1:
861 if x.size != y.size:
862 raise ValueError("Rankings must be of equal length.")
863 table = scipy.stats.contingency.crosstab(x, y)[1]
864 elif x.ndim == 2:
865 if np.any(x < 0):
866 raise ValueError("All elements of the contingency table must be "
867 "non-negative.")
868 if np.any(x != x.astype(int)):
869 raise ValueError("All elements of the contingency table must be "
870 "integer.")
871 if x.nonzero()[0].size < 2:
872 raise ValueError("At least two elements of the contingency table "
873 "must be nonzero.")
874 table = x
875 else:
876 raise ValueError("x must be either a 1D or 2D array")
877 # The table type is converted to a float to avoid an integer overflow
878 d, p = _somers_d(table.astype(float), alternative)
880 # add alias for consistency with other correlation functions
881 res = SomersDResult(d, p, table)
882 res.correlation = d
883 return res
886# This could be combined with `_all_partitions` in `_resampling.py`
887def _all_partitions(nx, ny):
888 """
889 Partition a set of indices into two fixed-length sets in all possible ways
891 Partition a set of indices 0 ... nx + ny - 1 into two sets of length nx and
892 ny in all possible ways (ignoring order of elements).
893 """
894 z = np.arange(nx+ny)
895 for c in combinations(z, nx):
896 x = np.array(c)
897 mask = np.ones(nx+ny, bool)
898 mask[x] = False
899 y = z[mask]
900 yield x, y
903def _compute_log_combinations(n):
904 """Compute all log combination of C(n, k)."""
905 gammaln_arr = gammaln(np.arange(n + 1) + 1)
906 return gammaln(n + 1) - gammaln_arr - gammaln_arr[::-1]
909@dataclass
910class BarnardExactResult:
911 statistic: float
912 pvalue: float
915def barnard_exact(table, alternative="two-sided", pooled=True, n=32):
916 r"""Perform a Barnard exact test on a 2x2 contingency table.
918 Parameters
919 ----------
920 table : array_like of ints
921 A 2x2 contingency table. Elements should be non-negative integers.
923 alternative : {'two-sided', 'less', 'greater'}, optional
924 Defines the null and alternative hypotheses. Default is 'two-sided'.
925 Please see explanations in the Notes section below.
927 pooled : bool, optional
928 Whether to compute score statistic with pooled variance (as in
929 Student's t-test, for example) or unpooled variance (as in Welch's
930 t-test). Default is ``True``.
932 n : int, optional
933 Number of sampling points used in the construction of the sampling
934 method. Note that this argument will automatically be converted to
935 the next higher power of 2 since `scipy.stats.qmc.Sobol` is used to
936 select sample points. Default is 32. Must be positive. In most cases,
937 32 points is enough to reach good precision. More points comes at
938 performance cost.
940 Returns
941 -------
942 ber : BarnardExactResult
943 A result object with the following attributes.
945 statistic : float
946 The Wald statistic with pooled or unpooled variance, depending
947 on the user choice of `pooled`.
949 pvalue : float
950 P-value, the probability of obtaining a distribution at least as
951 extreme as the one that was actually observed, assuming that the
952 null hypothesis is true.
954 See Also
955 --------
956 chi2_contingency : Chi-square test of independence of variables in a
957 contingency table.
958 fisher_exact : Fisher exact test on a 2x2 contingency table.
959 boschloo_exact : Boschloo's exact test on a 2x2 contingency table,
960 which is an uniformly more powerful alternative to Fisher's exact test.
962 Notes
963 -----
964 Barnard's test is an exact test used in the analysis of contingency
965 tables. It examines the association of two categorical variables, and
966 is a more powerful alternative than Fisher's exact test
967 for 2x2 contingency tables.
969 Let's define :math:`X_0` a 2x2 matrix representing the observed sample,
970 where each column stores the binomial experiment, as in the example
971 below. Let's also define :math:`p_1, p_2` the theoretical binomial
972 probabilities for :math:`x_{11}` and :math:`x_{12}`. When using
973 Barnard exact test, we can assert three different null hypotheses :
975 - :math:`H_0 : p_1 \geq p_2` versus :math:`H_1 : p_1 < p_2`,
976 with `alternative` = "less"
978 - :math:`H_0 : p_1 \leq p_2` versus :math:`H_1 : p_1 > p_2`,
979 with `alternative` = "greater"
981 - :math:`H_0 : p_1 = p_2` versus :math:`H_1 : p_1 \neq p_2`,
982 with `alternative` = "two-sided" (default one)
984 In order to compute Barnard's exact test, we are using the Wald
985 statistic [3]_ with pooled or unpooled variance.
986 Under the default assumption that both variances are equal
987 (``pooled = True``), the statistic is computed as:
989 .. math::
991 T(X) = \frac{
992 \hat{p}_1 - \hat{p}_2
993 }{
994 \sqrt{
995 \hat{p}(1 - \hat{p})
996 (\frac{1}{c_1} +
997 \frac{1}{c_2})
998 }
999 }
1001 with :math:`\hat{p}_1, \hat{p}_2` and :math:`\hat{p}` the estimator of
1002 :math:`p_1, p_2` and :math:`p`, the latter being the combined probability,
1003 given the assumption that :math:`p_1 = p_2`.
1005 If this assumption is invalid (``pooled = False``), the statistic is:
1007 .. math::
1009 T(X) = \frac{
1010 \hat{p}_1 - \hat{p}_2
1011 }{
1012 \sqrt{
1013 \frac{\hat{p}_1 (1 - \hat{p}_1)}{c_1} +
1014 \frac{\hat{p}_2 (1 - \hat{p}_2)}{c_2}
1015 }
1016 }
1018 The p-value is then computed as:
1020 .. math::
1022 \sum
1023 \binom{c_1}{x_{11}}
1024 \binom{c_2}{x_{12}}
1025 \pi^{x_{11} + x_{12}}
1026 (1 - \pi)^{t - x_{11} - x_{12}}
1028 where the sum is over all 2x2 contingency tables :math:`X` such that:
1029 * :math:`T(X) \leq T(X_0)` when `alternative` = "less",
1030 * :math:`T(X) \geq T(X_0)` when `alternative` = "greater", or
1031 * :math:`T(X) \geq |T(X_0)|` when `alternative` = "two-sided".
1032 Above, :math:`c_1, c_2` are the sum of the columns 1 and 2,
1033 and :math:`t` the total (sum of the 4 sample's element).
1035 The returned p-value is the maximum p-value taken over the nuisance
1036 parameter :math:`\pi`, where :math:`0 \leq \pi \leq 1`.
1038 This function's complexity is :math:`O(n c_1 c_2)`, where `n` is the
1039 number of sample points.
1041 References
1042 ----------
1043 .. [1] Barnard, G. A. "Significance Tests for 2x2 Tables". *Biometrika*.
1044 34.1/2 (1947): 123-138. :doi:`dpgkg3`
1046 .. [2] Mehta, Cyrus R., and Pralay Senchaudhuri. "Conditional versus
1047 unconditional exact tests for comparing two binomials."
1048 *Cytel Software Corporation* 675 (2003): 1-5.
1050 .. [3] "Wald Test". *Wikipedia*. https://en.wikipedia.org/wiki/Wald_test
1052 Examples
1053 --------
1054 An example use of Barnard's test is presented in [2]_.
1056 Consider the following example of a vaccine efficacy study
1057 (Chan, 1998). In a randomized clinical trial of 30 subjects, 15 were
1058 inoculated with a recombinant DNA influenza vaccine and the 15 were
1059 inoculated with a placebo. Twelve of the 15 subjects in the placebo
1060 group (80%) eventually became infected with influenza whereas for the
1061 vaccine group, only 7 of the 15 subjects (47%) became infected. The
1062 data are tabulated as a 2 x 2 table::
1064 Vaccine Placebo
1065 Yes 7 12
1066 No 8 3
1068 When working with statistical hypothesis testing, we usually use a
1069 threshold probability or significance level upon which we decide
1070 to reject the null hypothesis :math:`H_0`. Suppose we choose the common
1071 significance level of 5%.
1073 Our alternative hypothesis is that the vaccine will lower the chance of
1074 becoming infected with the virus; that is, the probability :math:`p_1` of
1075 catching the virus with the vaccine will be *less than* the probability
1076 :math:`p_2` of catching the virus without the vaccine. Therefore, we call
1077 `barnard_exact` with the ``alternative="less"`` option:
1079 >>> import scipy.stats as stats
1080 >>> res = stats.barnard_exact([[7, 12], [8, 3]], alternative="less")
1081 >>> res.statistic
1082 -1.894...
1083 >>> res.pvalue
1084 0.03407...
1086 Under the null hypothesis that the vaccine will not lower the chance of
1087 becoming infected, the probability of obtaining test results at least as
1088 extreme as the observed data is approximately 3.4%. Since this p-value is
1089 less than our chosen significance level, we have evidence to reject
1090 :math:`H_0` in favor of the alternative.
1092 Suppose we had used Fisher's exact test instead:
1094 >>> _, pvalue = stats.fisher_exact([[7, 12], [8, 3]], alternative="less")
1095 >>> pvalue
1096 0.0640...
1098 With the same threshold significance of 5%, we would not have been able
1099 to reject the null hypothesis in favor of the alternative. As stated in
1100 [2]_, Barnard's test is uniformly more powerful than Fisher's exact test
1101 because Barnard's test does not condition on any margin. Fisher's test
1102 should only be used when both sets of marginals are fixed.
1104 """
1105 if n <= 0:
1106 raise ValueError(
1107 "Number of points `n` must be strictly positive, "
1108 f"found {n!r}"
1109 )
1111 table = np.asarray(table, dtype=np.int64)
1113 if not table.shape == (2, 2):
1114 raise ValueError("The input `table` must be of shape (2, 2).")
1116 if np.any(table < 0):
1117 raise ValueError("All values in `table` must be nonnegative.")
1119 if 0 in table.sum(axis=0):
1120 # If both values in column are zero, the p-value is 1 and
1121 # the score's statistic is NaN.
1122 return BarnardExactResult(np.nan, 1.0)
1124 total_col_1, total_col_2 = table.sum(axis=0)
1126 x1 = np.arange(total_col_1 + 1, dtype=np.int64).reshape(-1, 1)
1127 x2 = np.arange(total_col_2 + 1, dtype=np.int64).reshape(1, -1)
1129 # We need to calculate the wald statistics for each combination of x1 and
1130 # x2.
1131 p1, p2 = x1 / total_col_1, x2 / total_col_2
1133 if pooled:
1134 p = (x1 + x2) / (total_col_1 + total_col_2)
1135 variances = p * (1 - p) * (1 / total_col_1 + 1 / total_col_2)
1136 else:
1137 variances = p1 * (1 - p1) / total_col_1 + p2 * (1 - p2) / total_col_2
1139 # To avoid warning when dividing by 0
1140 with np.errstate(divide="ignore", invalid="ignore"):
1141 wald_statistic = np.divide((p1 - p2), np.sqrt(variances))
1143 wald_statistic[p1 == p2] = 0 # Removing NaN values
1145 wald_stat_obs = wald_statistic[table[0, 0], table[0, 1]]
1147 if alternative == "two-sided":
1148 index_arr = np.abs(wald_statistic) >= abs(wald_stat_obs)
1149 elif alternative == "less":
1150 index_arr = wald_statistic <= wald_stat_obs
1151 elif alternative == "greater":
1152 index_arr = wald_statistic >= wald_stat_obs
1153 else:
1154 msg = (
1155 "`alternative` should be one of {'two-sided', 'less', 'greater'},"
1156 f" found {alternative!r}"
1157 )
1158 raise ValueError(msg)
1160 x1_sum_x2 = x1 + x2
1162 x1_log_comb = _compute_log_combinations(total_col_1)
1163 x2_log_comb = _compute_log_combinations(total_col_2)
1164 x1_sum_x2_log_comb = x1_log_comb[x1] + x2_log_comb[x2]
1166 result = shgo(
1167 _get_binomial_log_p_value_with_nuisance_param,
1168 args=(x1_sum_x2, x1_sum_x2_log_comb, index_arr),
1169 bounds=((0, 1),),
1170 n=n,
1171 sampling_method="sobol",
1172 )
1174 # result.fun is the negative log pvalue and therefore needs to be
1175 # changed before return
1176 p_value = np.clip(np.exp(-result.fun), a_min=0, a_max=1)
1177 return BarnardExactResult(wald_stat_obs, p_value)
1180@dataclass
1181class BoschlooExactResult:
1182 statistic: float
1183 pvalue: float
1186def boschloo_exact(table, alternative="two-sided", n=32):
1187 r"""Perform Boschloo's exact test on a 2x2 contingency table.
1189 Parameters
1190 ----------
1191 table : array_like of ints
1192 A 2x2 contingency table. Elements should be non-negative integers.
1194 alternative : {'two-sided', 'less', 'greater'}, optional
1195 Defines the null and alternative hypotheses. Default is 'two-sided'.
1196 Please see explanations in the Notes section below.
1198 n : int, optional
1199 Number of sampling points used in the construction of the sampling
1200 method. Note that this argument will automatically be converted to
1201 the next higher power of 2 since `scipy.stats.qmc.Sobol` is used to
1202 select sample points. Default is 32. Must be positive. In most cases,
1203 32 points is enough to reach good precision. More points comes at
1204 performance cost.
1206 Returns
1207 -------
1208 ber : BoschlooExactResult
1209 A result object with the following attributes.
1211 statistic : float
1212 The statistic used in Boschloo's test; that is, the p-value
1213 from Fisher's exact test.
1215 pvalue : float
1216 P-value, the probability of obtaining a distribution at least as
1217 extreme as the one that was actually observed, assuming that the
1218 null hypothesis is true.
1220 See Also
1221 --------
1222 chi2_contingency : Chi-square test of independence of variables in a
1223 contingency table.
1224 fisher_exact : Fisher exact test on a 2x2 contingency table.
1225 barnard_exact : Barnard's exact test, which is a more powerful alternative
1226 than Fisher's exact test for 2x2 contingency tables.
1228 Notes
1229 -----
1230 Boschloo's test is an exact test used in the analysis of contingency
1231 tables. It examines the association of two categorical variables, and
1232 is a uniformly more powerful alternative to Fisher's exact test
1233 for 2x2 contingency tables.
1235 Boschloo's exact test uses the p-value of Fisher's exact test as a
1236 statistic, and Boschloo's p-value is the probability under the null
1237 hypothesis of observing such an extreme value of this statistic.
1239 Let's define :math:`X_0` a 2x2 matrix representing the observed sample,
1240 where each column stores the binomial experiment, as in the example
1241 below. Let's also define :math:`p_1, p_2` the theoretical binomial
1242 probabilities for :math:`x_{11}` and :math:`x_{12}`. When using
1243 Boschloo exact test, we can assert three different alternative hypotheses:
1245 - :math:`H_0 : p_1=p_2` versus :math:`H_1 : p_1 < p_2`,
1246 with `alternative` = "less"
1248 - :math:`H_0 : p_1=p_2` versus :math:`H_1 : p_1 > p_2`,
1249 with `alternative` = "greater"
1251 - :math:`H_0 : p_1=p_2` versus :math:`H_1 : p_1 \neq p_2`,
1252 with `alternative` = "two-sided" (default)
1254 There are multiple conventions for computing a two-sided p-value when the
1255 null distribution is asymmetric. Here, we apply the convention that the
1256 p-value of a two-sided test is twice the minimum of the p-values of the
1257 one-sided tests (clipped to 1.0). Note that `fisher_exact` follows a
1258 different convention, so for a given `table`, the statistic reported by
1259 `boschloo_exact` may differ from the p-value reported by `fisher_exact`
1260 when ``alternative='two-sided'``.
1262 .. versionadded:: 1.7.0
1264 References
1265 ----------
1266 .. [1] R.D. Boschloo. "Raised conditional level of significance for the
1267 2 x 2-table when testing the equality of two probabilities",
1268 Statistica Neerlandica, 24(1), 1970
1270 .. [2] "Boschloo's test", Wikipedia,
1271 https://en.wikipedia.org/wiki/Boschloo%27s_test
1273 .. [3] Lise M. Saari et al. "Employee attitudes and job satisfaction",
1274 Human Resource Management, 43(4), 395-407, 2004,
1275 :doi:`10.1002/hrm.20032`.
1277 Examples
1278 --------
1279 In the following example, we consider the article "Employee
1280 attitudes and job satisfaction" [3]_
1281 which reports the results of a survey from 63 scientists and 117 college
1282 professors. Of the 63 scientists, 31 said they were very satisfied with
1283 their jobs, whereas 74 of the college professors were very satisfied
1284 with their work. Is this significant evidence that college
1285 professors are happier with their work than scientists?
1286 The following table summarizes the data mentioned above::
1288 college professors scientists
1289 Very Satisfied 74 31
1290 Dissatisfied 43 32
1292 When working with statistical hypothesis testing, we usually use a
1293 threshold probability or significance level upon which we decide
1294 to reject the null hypothesis :math:`H_0`. Suppose we choose the common
1295 significance level of 5%.
1297 Our alternative hypothesis is that college professors are truly more
1298 satisfied with their work than scientists. Therefore, we expect
1299 :math:`p_1` the proportion of very satisfied college professors to be
1300 greater than :math:`p_2`, the proportion of very satisfied scientists.
1301 We thus call `boschloo_exact` with the ``alternative="greater"`` option:
1303 >>> import scipy.stats as stats
1304 >>> res = stats.boschloo_exact([[74, 31], [43, 32]], alternative="greater")
1305 >>> res.statistic
1306 0.0483...
1307 >>> res.pvalue
1308 0.0355...
1310 Under the null hypothesis that scientists are happier in their work than
1311 college professors, the probability of obtaining test
1312 results at least as extreme as the observed data is approximately 3.55%.
1313 Since this p-value is less than our chosen significance level, we have
1314 evidence to reject :math:`H_0` in favor of the alternative hypothesis.
1316 """
1317 hypergeom = distributions.hypergeom
1319 if n <= 0:
1320 raise ValueError(
1321 "Number of points `n` must be strictly positive,"
1322 f" found {n!r}"
1323 )
1325 table = np.asarray(table, dtype=np.int64)
1327 if not table.shape == (2, 2):
1328 raise ValueError("The input `table` must be of shape (2, 2).")
1330 if np.any(table < 0):
1331 raise ValueError("All values in `table` must be nonnegative.")
1333 if 0 in table.sum(axis=0):
1334 # If both values in column are zero, the p-value is 1 and
1335 # the score's statistic is NaN.
1336 return BoschlooExactResult(np.nan, np.nan)
1338 total_col_1, total_col_2 = table.sum(axis=0)
1339 total = total_col_1 + total_col_2
1340 x1 = np.arange(total_col_1 + 1, dtype=np.int64).reshape(1, -1)
1341 x2 = np.arange(total_col_2 + 1, dtype=np.int64).reshape(-1, 1)
1342 x1_sum_x2 = x1 + x2
1344 if alternative == 'less':
1345 pvalues = hypergeom.cdf(x1, total, x1_sum_x2, total_col_1).T
1346 elif alternative == 'greater':
1347 # Same formula as the 'less' case, but with the second column.
1348 pvalues = hypergeom.cdf(x2, total, x1_sum_x2, total_col_2).T
1349 elif alternative == 'two-sided':
1350 boschloo_less = boschloo_exact(table, alternative="less", n=n)
1351 boschloo_greater = boschloo_exact(table, alternative="greater", n=n)
1353 res = (
1354 boschloo_less if boschloo_less.pvalue < boschloo_greater.pvalue
1355 else boschloo_greater
1356 )
1358 # Two-sided p-value is defined as twice the minimum of the one-sided
1359 # p-values
1360 pvalue = np.clip(2 * res.pvalue, a_min=0, a_max=1)
1361 return BoschlooExactResult(res.statistic, pvalue)
1362 else:
1363 msg = (
1364 f"`alternative` should be one of {'two-sided', 'less', 'greater'},"
1365 f" found {alternative!r}"
1366 )
1367 raise ValueError(msg)
1369 fisher_stat = pvalues[table[0, 0], table[0, 1]]
1371 # fisher_stat * (1+1e-13) guards us from small numerical error. It is
1372 # equivalent to np.isclose with relative tol of 1e-13 and absolute tol of 0
1373 # For more throughout explanations, see gh-14178
1374 index_arr = pvalues <= fisher_stat * (1+1e-13)
1376 x1, x2, x1_sum_x2 = x1.T, x2.T, x1_sum_x2.T
1377 x1_log_comb = _compute_log_combinations(total_col_1)
1378 x2_log_comb = _compute_log_combinations(total_col_2)
1379 x1_sum_x2_log_comb = x1_log_comb[x1] + x2_log_comb[x2]
1381 result = shgo(
1382 _get_binomial_log_p_value_with_nuisance_param,
1383 args=(x1_sum_x2, x1_sum_x2_log_comb, index_arr),
1384 bounds=((0, 1),),
1385 n=n,
1386 sampling_method="sobol",
1387 )
1389 # result.fun is the negative log pvalue and therefore needs to be
1390 # changed before return
1391 p_value = np.clip(np.exp(-result.fun), a_min=0, a_max=1)
1392 return BoschlooExactResult(fisher_stat, p_value)
1395def _get_binomial_log_p_value_with_nuisance_param(
1396 nuisance_param, x1_sum_x2, x1_sum_x2_log_comb, index_arr
1397):
1398 r"""
1399 Compute the log pvalue in respect of a nuisance parameter considering
1400 a 2x2 sample space.
1402 Parameters
1403 ----------
1404 nuisance_param : float
1405 nuisance parameter used in the computation of the maximisation of
1406 the p-value. Must be between 0 and 1
1408 x1_sum_x2 : ndarray
1409 Sum of x1 and x2 inside barnard_exact
1411 x1_sum_x2_log_comb : ndarray
1412 sum of the log combination of x1 and x2
1414 index_arr : ndarray of boolean
1416 Returns
1417 -------
1418 p_value : float
1419 Return the maximum p-value considering every nuisance paramater
1420 between 0 and 1
1422 Notes
1423 -----
1425 Both Barnard's test and Boschloo's test iterate over a nuisance parameter
1426 :math:`\pi \in [0, 1]` to find the maximum p-value. To search this
1427 maxima, this function return the negative log pvalue with respect to the
1428 nuisance parameter passed in params. This negative log p-value is then
1429 used in `shgo` to find the minimum negative pvalue which is our maximum
1430 pvalue.
1432 Also, to compute the different combination used in the
1433 p-values' computation formula, this function uses `gammaln` which is
1434 more tolerant for large value than `scipy.special.comb`. `gammaln` gives
1435 a log combination. For the little precision loss, performances are
1436 improved a lot.
1437 """
1438 t1, t2 = x1_sum_x2.shape
1439 n = t1 + t2 - 2
1440 with np.errstate(divide="ignore", invalid="ignore"):
1441 log_nuisance = np.log(
1442 nuisance_param,
1443 out=np.zeros_like(nuisance_param),
1444 where=nuisance_param >= 0,
1445 )
1446 log_1_minus_nuisance = np.log(
1447 1 - nuisance_param,
1448 out=np.zeros_like(nuisance_param),
1449 where=1 - nuisance_param >= 0,
1450 )
1452 nuisance_power_x1_x2 = log_nuisance * x1_sum_x2
1453 nuisance_power_x1_x2[(x1_sum_x2 == 0)[:, :]] = 0
1455 nuisance_power_n_minus_x1_x2 = log_1_minus_nuisance * (n - x1_sum_x2)
1456 nuisance_power_n_minus_x1_x2[(x1_sum_x2 == n)[:, :]] = 0
1458 tmp_log_values_arr = (
1459 x1_sum_x2_log_comb
1460 + nuisance_power_x1_x2
1461 + nuisance_power_n_minus_x1_x2
1462 )
1464 tmp_values_from_index = tmp_log_values_arr[index_arr]
1466 # To avoid dividing by zero in log function and getting inf value,
1467 # values are centered according to the max
1468 max_value = tmp_values_from_index.max()
1470 # To have better result's precision, the log pvalue is taken here.
1471 # Indeed, pvalue is included inside [0, 1] interval. Passing the
1472 # pvalue to log makes the interval a lot bigger ([-inf, 0]), and thus
1473 # help us to achieve better precision
1474 with np.errstate(divide="ignore", invalid="ignore"):
1475 log_probs = np.exp(tmp_values_from_index - max_value).sum()
1476 log_pvalue = max_value + np.log(
1477 log_probs,
1478 out=np.full_like(log_probs, -np.inf),
1479 where=log_probs > 0,
1480 )
1482 # Since shgo find the minima, minus log pvalue is returned
1483 return -log_pvalue
1486def _pval_cvm_2samp_exact(s, m, n):
1487 """
1488 Compute the exact p-value of the Cramer-von Mises two-sample test
1489 for a given value s of the test statistic.
1490 m and n are the sizes of the samples.
1492 [1] Y. Xiao, A. Gordon, and A. Yakovlev, "A C++ Program for
1493 the Cramér-Von Mises Two-Sample Test", J. Stat. Soft.,
1494 vol. 17, no. 8, pp. 1-15, Dec. 2006.
1495 [2] T. W. Anderson "On the Distribution of the Two-Sample Cramer-von Mises
1496 Criterion," The Annals of Mathematical Statistics, Ann. Math. Statist.
1497 33(3), 1148-1159, (September, 1962)
1498 """
1500 # [1, p. 3]
1501 lcm = np.lcm(m, n)
1502 # [1, p. 4], below eq. 3
1503 a = lcm // m
1504 b = lcm // n
1505 # Combine Eq. 9 in [2] with Eq. 2 in [1] and solve for $\zeta$
1506 # Hint: `s` is $U$ in [2], and $T_2$ in [1] is $T$ in [2]
1507 mn = m * n
1508 zeta = lcm ** 2 * (m + n) * (6 * s - mn * (4 * mn - 1)) // (6 * mn ** 2)
1510 # bound maximum value that may appear in `gs` (remember both rows!)
1511 zeta_bound = lcm**2 * (m + n) # bound elements in row 1
1512 combinations = comb(m + n, m) # sum of row 2
1513 max_gs = max(zeta_bound, combinations)
1514 dtype = np.min_scalar_type(max_gs)
1516 # the frequency table of $g_{u, v}^+$ defined in [1, p. 6]
1517 gs = ([np.array([[0], [1]], dtype=dtype)]
1518 + [np.empty((2, 0), dtype=dtype) for _ in range(m)])
1519 for u in range(n + 1):
1520 next_gs = []
1521 tmp = np.empty((2, 0), dtype=dtype)
1522 for v, g in enumerate(gs):
1523 # Calculate g recursively with eq. 11 in [1]. Even though it
1524 # doesn't look like it, this also does 12/13 (all of Algorithm 1).
1525 vi, i0, i1 = np.intersect1d(tmp[0], g[0], return_indices=True)
1526 tmp = np.concatenate([
1527 np.stack([vi, tmp[1, i0] + g[1, i1]]),
1528 np.delete(tmp, i0, 1),
1529 np.delete(g, i1, 1)
1530 ], 1)
1531 tmp[0] += (a * v - b * u) ** 2
1532 next_gs.append(tmp)
1533 gs = next_gs
1534 value, freq = gs[m]
1535 return np.float64(np.sum(freq[value >= zeta]) / combinations)
1538def cramervonmises_2samp(x, y, method='auto'):
1539 """Perform the two-sample Cramér-von Mises test for goodness of fit.
1541 This is the two-sample version of the Cramér-von Mises test ([1]_):
1542 for two independent samples :math:`X_1, ..., X_n` and
1543 :math:`Y_1, ..., Y_m`, the null hypothesis is that the samples
1544 come from the same (unspecified) continuous distribution.
1546 Parameters
1547 ----------
1548 x : array_like
1549 A 1-D array of observed values of the random variables :math:`X_i`.
1550 y : array_like
1551 A 1-D array of observed values of the random variables :math:`Y_i`.
1552 method : {'auto', 'asymptotic', 'exact'}, optional
1553 The method used to compute the p-value, see Notes for details.
1554 The default is 'auto'.
1556 Returns
1557 -------
1558 res : object with attributes
1559 statistic : float
1560 Cramér-von Mises statistic.
1561 pvalue : float
1562 The p-value.
1564 See Also
1565 --------
1566 cramervonmises, anderson_ksamp, epps_singleton_2samp, ks_2samp
1568 Notes
1569 -----
1570 .. versionadded:: 1.7.0
1572 The statistic is computed according to equation 9 in [2]_. The
1573 calculation of the p-value depends on the keyword `method`:
1575 - ``asymptotic``: The p-value is approximated by using the limiting
1576 distribution of the test statistic.
1577 - ``exact``: The exact p-value is computed by enumerating all
1578 possible combinations of the test statistic, see [2]_.
1580 If ``method='auto'``, the exact approach is used
1581 if both samples contain equal to or less than 20 observations,
1582 otherwise the asymptotic distribution is used.
1584 If the underlying distribution is not continuous, the p-value is likely to
1585 be conservative (Section 6.2 in [3]_). When ranking the data to compute
1586 the test statistic, midranks are used if there are ties.
1588 References
1589 ----------
1590 .. [1] https://en.wikipedia.org/wiki/Cramer-von_Mises_criterion
1591 .. [2] Anderson, T.W. (1962). On the distribution of the two-sample
1592 Cramer-von-Mises criterion. The Annals of Mathematical
1593 Statistics, pp. 1148-1159.
1594 .. [3] Conover, W.J., Practical Nonparametric Statistics, 1971.
1596 Examples
1597 --------
1599 Suppose we wish to test whether two samples generated by
1600 ``scipy.stats.norm.rvs`` have the same distribution. We choose a
1601 significance level of alpha=0.05.
1603 >>> import numpy as np
1604 >>> from scipy import stats
1605 >>> rng = np.random.default_rng()
1606 >>> x = stats.norm.rvs(size=100, random_state=rng)
1607 >>> y = stats.norm.rvs(size=70, random_state=rng)
1608 >>> res = stats.cramervonmises_2samp(x, y)
1609 >>> res.statistic, res.pvalue
1610 (0.29376470588235293, 0.1412873014573014)
1612 The p-value exceeds our chosen significance level, so we do not
1613 reject the null hypothesis that the observed samples are drawn from the
1614 same distribution.
1616 For small sample sizes, one can compute the exact p-values:
1618 >>> x = stats.norm.rvs(size=7, random_state=rng)
1619 >>> y = stats.t.rvs(df=2, size=6, random_state=rng)
1620 >>> res = stats.cramervonmises_2samp(x, y, method='exact')
1621 >>> res.statistic, res.pvalue
1622 (0.197802197802198, 0.31643356643356646)
1624 The p-value based on the asymptotic distribution is a good approximation
1625 even though the sample size is small.
1627 >>> res = stats.cramervonmises_2samp(x, y, method='asymptotic')
1628 >>> res.statistic, res.pvalue
1629 (0.197802197802198, 0.2966041181527128)
1631 Independent of the method, one would not reject the null hypothesis at the
1632 chosen significance level in this example.
1634 """
1635 xa = np.sort(np.asarray(x))
1636 ya = np.sort(np.asarray(y))
1638 if xa.size <= 1 or ya.size <= 1:
1639 raise ValueError('x and y must contain at least two observations.')
1640 if xa.ndim > 1 or ya.ndim > 1:
1641 raise ValueError('The samples must be one-dimensional.')
1642 if method not in ['auto', 'exact', 'asymptotic']:
1643 raise ValueError('method must be either auto, exact or asymptotic.')
1645 nx = len(xa)
1646 ny = len(ya)
1648 if method == 'auto':
1649 if max(nx, ny) > 20:
1650 method = 'asymptotic'
1651 else:
1652 method = 'exact'
1654 # get ranks of x and y in the pooled sample
1655 z = np.concatenate([xa, ya])
1656 # in case of ties, use midrank (see [1])
1657 r = scipy.stats.rankdata(z, method='average')
1658 rx = r[:nx]
1659 ry = r[nx:]
1661 # compute U (eq. 10 in [2])
1662 u = nx * np.sum((rx - np.arange(1, nx+1))**2)
1663 u += ny * np.sum((ry - np.arange(1, ny+1))**2)
1665 # compute T (eq. 9 in [2])
1666 k, N = nx*ny, nx + ny
1667 t = u / (k*N) - (4*k - 1)/(6*N)
1669 if method == 'exact':
1670 p = _pval_cvm_2samp_exact(u, nx, ny)
1671 else:
1672 # compute expected value and variance of T (eq. 11 and 14 in [2])
1673 et = (1 + 1/N)/6
1674 vt = (N+1) * (4*k*N - 3*(nx**2 + ny**2) - 2*k)
1675 vt = vt / (45 * N**2 * 4 * k)
1677 # computed the normalized statistic (eq. 15 in [2])
1678 tn = 1/6 + (t - et) / np.sqrt(45 * vt)
1680 # approximate distribution of tn with limiting distribution
1681 # of the one-sample test statistic
1682 # if tn < 0.003, the _cdf_cvm_inf(tn) < 1.28*1e-18, return 1.0 directly
1683 if tn < 0.003:
1684 p = 1.0
1685 else:
1686 p = max(0, 1. - _cdf_cvm_inf(tn))
1688 return CramerVonMisesResult(statistic=t, pvalue=p)
1691class TukeyHSDResult:
1692 """Result of `scipy.stats.tukey_hsd`.
1694 Attributes
1695 ----------
1696 statistic : float ndarray
1697 The computed statistic of the test for each comparison. The element
1698 at index ``(i, j)`` is the statistic for the comparison between groups
1699 ``i`` and ``j``.
1700 pvalue : float ndarray
1701 The associated p-value from the studentized range distribution. The
1702 element at index ``(i, j)`` is the p-value for the comparison
1703 between groups ``i`` and ``j``.
1705 Notes
1706 -----
1707 The string representation of this object displays the most recently
1708 calculated confidence interval, and if none have been previously
1709 calculated, it will evaluate ``confidence_interval()``.
1711 References
1712 ----------
1713 .. [1] NIST/SEMATECH e-Handbook of Statistical Methods, "7.4.7.1. Tukey's
1714 Method."
1715 https://www.itl.nist.gov/div898/handbook/prc/section4/prc471.htm,
1716 28 November 2020.
1717 """
1719 def __init__(self, statistic, pvalue, _nobs, _ntreatments, _stand_err):
1720 self.statistic = statistic
1721 self.pvalue = pvalue
1722 self._ntreatments = _ntreatments
1723 self._nobs = _nobs
1724 self._stand_err = _stand_err
1725 self._ci = None
1726 self._ci_cl = None
1728 def __str__(self):
1729 # Note: `__str__` prints the confidence intervals from the most
1730 # recent call to `confidence_interval`. If it has not been called,
1731 # it will be called with the default CL of .95.
1732 if self._ci is None:
1733 self.confidence_interval(confidence_level=.95)
1734 s = ("Tukey's HSD Pairwise Group Comparisons"
1735 f" ({self._ci_cl*100:.1f}% Confidence Interval)\n")
1736 s += "Comparison Statistic p-value Lower CI Upper CI\n"
1737 for i in range(self.pvalue.shape[0]):
1738 for j in range(self.pvalue.shape[0]):
1739 if i != j:
1740 s += (f" ({i} - {j}) {self.statistic[i, j]:>10.3f}"
1741 f"{self.pvalue[i, j]:>10.3f}"
1742 f"{self._ci.low[i, j]:>10.3f}"
1743 f"{self._ci.high[i, j]:>10.3f}\n")
1744 return s
1746 def confidence_interval(self, confidence_level=.95):
1747 """Compute the confidence interval for the specified confidence level.
1749 Parameters
1750 ----------
1751 confidence_level : float, optional
1752 Confidence level for the computed confidence interval
1753 of the estimated proportion. Default is .95.
1755 Returns
1756 -------
1757 ci : ``ConfidenceInterval`` object
1758 The object has attributes ``low`` and ``high`` that hold the
1759 lower and upper bounds of the confidence intervals for each
1760 comparison. The high and low values are accessible for each
1761 comparison at index ``(i, j)`` between groups ``i`` and ``j``.
1763 References
1764 ----------
1765 .. [1] NIST/SEMATECH e-Handbook of Statistical Methods, "7.4.7.1.
1766 Tukey's Method."
1767 https://www.itl.nist.gov/div898/handbook/prc/section4/prc471.htm,
1768 28 November 2020.
1770 Examples
1771 --------
1772 >>> from scipy.stats import tukey_hsd
1773 >>> group0 = [24.5, 23.5, 26.4, 27.1, 29.9]
1774 >>> group1 = [28.4, 34.2, 29.5, 32.2, 30.1]
1775 >>> group2 = [26.1, 28.3, 24.3, 26.2, 27.8]
1776 >>> result = tukey_hsd(group0, group1, group2)
1777 >>> ci = result.confidence_interval()
1778 >>> ci.low
1779 array([[-3.649159, -8.249159, -3.909159],
1780 [ 0.950841, -3.649159, 0.690841],
1781 [-3.389159, -7.989159, -3.649159]])
1782 >>> ci.high
1783 array([[ 3.649159, -0.950841, 3.389159],
1784 [ 8.249159, 3.649159, 7.989159],
1785 [ 3.909159, -0.690841, 3.649159]])
1786 """
1787 # check to see if the supplied confidence level matches that of the
1788 # previously computed CI.
1789 if (self._ci is not None and self._ci_cl is not None and
1790 confidence_level == self._ci_cl):
1791 return self._ci
1793 if not 0 < confidence_level < 1:
1794 raise ValueError("Confidence level must be between 0 and 1.")
1795 # determine the critical value of the studentized range using the
1796 # appropriate confidence level, number of treatments, and degrees
1797 # of freedom as determined by the number of data less the number of
1798 # treatments. ("Confidence limits for Tukey's method")[1]. Note that
1799 # in the cases of unequal sample sizes there will be a criterion for
1800 # each group comparison.
1801 params = (confidence_level, self._nobs, self._ntreatments - self._nobs)
1802 srd = distributions.studentized_range.ppf(*params)
1803 # also called maximum critical value, the Tukey criterion is the
1804 # studentized range critical value * the square root of mean square
1805 # error over the sample size.
1806 tukey_criterion = srd * self._stand_err
1807 # the confidence levels are determined by the
1808 # `mean_differences` +- `tukey_criterion`
1809 upper_conf = self.statistic + tukey_criterion
1810 lower_conf = self.statistic - tukey_criterion
1811 self._ci = ConfidenceInterval(low=lower_conf, high=upper_conf)
1812 self._ci_cl = confidence_level
1813 return self._ci
1816def _tukey_hsd_iv(args):
1817 if (len(args)) < 2:
1818 raise ValueError("There must be more than 1 treatment.")
1819 args = [np.asarray(arg) for arg in args]
1820 for arg in args:
1821 if arg.ndim != 1:
1822 raise ValueError("Input samples must be one-dimensional.")
1823 if arg.size <= 1:
1824 raise ValueError("Input sample size must be greater than one.")
1825 if np.isinf(arg).any():
1826 raise ValueError("Input samples must be finite.")
1827 return args
1830def tukey_hsd(*args):
1831 """Perform Tukey's HSD test for equality of means over multiple treatments.
1833 Tukey's honestly significant difference (HSD) test performs pairwise
1834 comparison of means for a set of samples. Whereas ANOVA (e.g. `f_oneway`)
1835 assesses whether the true means underlying each sample are identical,
1836 Tukey's HSD is a post hoc test used to compare the mean of each sample
1837 to the mean of each other sample.
1839 The null hypothesis is that the distributions underlying the samples all
1840 have the same mean. The test statistic, which is computed for every
1841 possible pairing of samples, is simply the difference between the sample
1842 means. For each pair, the p-value is the probability under the null
1843 hypothesis (and other assumptions; see notes) of observing such an extreme
1844 value of the statistic, considering that many pairwise comparisons are
1845 being performed. Confidence intervals for the difference between each pair
1846 of means are also available.
1848 Parameters
1849 ----------
1850 sample1, sample2, ... : array_like
1851 The sample measurements for each group. There must be at least
1852 two arguments.
1854 Returns
1855 -------
1856 result : `~scipy.stats._result_classes.TukeyHSDResult` instance
1857 The return value is an object with the following attributes:
1859 statistic : float ndarray
1860 The computed statistic of the test for each comparison. The element
1861 at index ``(i, j)`` is the statistic for the comparison between
1862 groups ``i`` and ``j``.
1863 pvalue : float ndarray
1864 The computed p-value of the test for each comparison. The element
1865 at index ``(i, j)`` is the p-value for the comparison between
1866 groups ``i`` and ``j``.
1868 The object has the following methods:
1870 confidence_interval(confidence_level=0.95):
1871 Compute the confidence interval for the specified confidence level.
1873 See Also
1874 --------
1875 dunnett : performs comparison of means against a control group.
1877 Notes
1878 -----
1879 The use of this test relies on several assumptions.
1881 1. The observations are independent within and among groups.
1882 2. The observations within each group are normally distributed.
1883 3. The distributions from which the samples are drawn have the same finite
1884 variance.
1886 The original formulation of the test was for samples of equal size [6]_.
1887 In case of unequal sample sizes, the test uses the Tukey-Kramer method
1888 [4]_.
1890 References
1891 ----------
1892 .. [1] NIST/SEMATECH e-Handbook of Statistical Methods, "7.4.7.1. Tukey's
1893 Method."
1894 https://www.itl.nist.gov/div898/handbook/prc/section4/prc471.htm,
1895 28 November 2020.
1896 .. [2] Abdi, Herve & Williams, Lynne. (2021). "Tukey's Honestly Significant
1897 Difference (HSD) Test."
1898 https://personal.utdallas.edu/~herve/abdi-HSD2010-pretty.pdf
1899 .. [3] "One-Way ANOVA Using SAS PROC ANOVA & PROC GLM." SAS
1900 Tutorials, 2007, www.stattutorials.com/SAS/TUTORIAL-PROC-GLM.htm.
1901 .. [4] Kramer, Clyde Young. "Extension of Multiple Range Tests to Group
1902 Means with Unequal Numbers of Replications." Biometrics, vol. 12,
1903 no. 3, 1956, pp. 307-310. JSTOR, www.jstor.org/stable/3001469.
1904 Accessed 25 May 2021.
1905 .. [5] NIST/SEMATECH e-Handbook of Statistical Methods, "7.4.3.3.
1906 The ANOVA table and tests of hypotheses about means"
1907 https://www.itl.nist.gov/div898/handbook/prc/section4/prc433.htm,
1908 2 June 2021.
1909 .. [6] Tukey, John W. "Comparing Individual Means in the Analysis of
1910 Variance." Biometrics, vol. 5, no. 2, 1949, pp. 99-114. JSTOR,
1911 www.jstor.org/stable/3001913. Accessed 14 June 2021.
1914 Examples
1915 --------
1916 Here are some data comparing the time to relief of three brands of
1917 headache medicine, reported in minutes. Data adapted from [3]_.
1919 >>> import numpy as np
1920 >>> from scipy.stats import tukey_hsd
1921 >>> group0 = [24.5, 23.5, 26.4, 27.1, 29.9]
1922 >>> group1 = [28.4, 34.2, 29.5, 32.2, 30.1]
1923 >>> group2 = [26.1, 28.3, 24.3, 26.2, 27.8]
1925 We would like to see if the means between any of the groups are
1926 significantly different. First, visually examine a box and whisker plot.
1928 >>> import matplotlib.pyplot as plt
1929 >>> fig, ax = plt.subplots(1, 1)
1930 >>> ax.boxplot([group0, group1, group2])
1931 >>> ax.set_xticklabels(["group0", "group1", "group2"]) # doctest: +SKIP
1932 >>> ax.set_ylabel("mean") # doctest: +SKIP
1933 >>> plt.show()
1935 From the box and whisker plot, we can see overlap in the interquartile
1936 ranges group 1 to group 2 and group 3, but we can apply the ``tukey_hsd``
1937 test to determine if the difference between means is significant. We
1938 set a significance level of .05 to reject the null hypothesis.
1940 >>> res = tukey_hsd(group0, group1, group2)
1941 >>> print(res)
1942 Tukey's HSD Pairwise Group Comparisons (95.0% Confidence Interval)
1943 Comparison Statistic p-value Lower CI Upper CI
1944 (0 - 1) -4.600 0.014 -8.249 -0.951
1945 (0 - 2) -0.260 0.980 -3.909 3.389
1946 (1 - 0) 4.600 0.014 0.951 8.249
1947 (1 - 2) 4.340 0.020 0.691 7.989
1948 (2 - 0) 0.260 0.980 -3.389 3.909
1949 (2 - 1) -4.340 0.020 -7.989 -0.691
1951 The null hypothesis is that each group has the same mean. The p-value for
1952 comparisons between ``group0`` and ``group1`` as well as ``group1`` and
1953 ``group2`` do not exceed .05, so we reject the null hypothesis that they
1954 have the same means. The p-value of the comparison between ``group0``
1955 and ``group2`` exceeds .05, so we accept the null hypothesis that there
1956 is not a significant difference between their means.
1958 We can also compute the confidence interval associated with our chosen
1959 confidence level.
1961 >>> group0 = [24.5, 23.5, 26.4, 27.1, 29.9]
1962 >>> group1 = [28.4, 34.2, 29.5, 32.2, 30.1]
1963 >>> group2 = [26.1, 28.3, 24.3, 26.2, 27.8]
1964 >>> result = tukey_hsd(group0, group1, group2)
1965 >>> conf = res.confidence_interval(confidence_level=.99)
1966 >>> for ((i, j), l) in np.ndenumerate(conf.low):
1967 ... # filter out self comparisons
1968 ... if i != j:
1969 ... h = conf.high[i,j]
1970 ... print(f"({i} - {j}) {l:>6.3f} {h:>6.3f}")
1971 (0 - 1) -9.480 0.280
1972 (0 - 2) -5.140 4.620
1973 (1 - 0) -0.280 9.480
1974 (1 - 2) -0.540 9.220
1975 (2 - 0) -4.620 5.140
1976 (2 - 1) -9.220 0.540
1977 """
1978 args = _tukey_hsd_iv(args)
1979 ntreatments = len(args)
1980 means = np.asarray([np.mean(arg) for arg in args])
1981 nsamples_treatments = np.asarray([a.size for a in args])
1982 nobs = np.sum(nsamples_treatments)
1984 # determine mean square error [5]. Note that this is sometimes called
1985 # mean square error within.
1986 mse = (np.sum([np.var(arg, ddof=1) for arg in args] *
1987 (nsamples_treatments - 1)) / (nobs - ntreatments))
1989 # The calculation of the standard error differs when treatments differ in
1990 # size. See ("Unequal sample sizes")[1].
1991 if np.unique(nsamples_treatments).size == 1:
1992 # all input groups are the same length, so only one value needs to be
1993 # calculated [1].
1994 normalize = 2 / nsamples_treatments[0]
1995 else:
1996 # to compare groups of differing sizes, we must compute a variance
1997 # value for each individual comparison. Use broadcasting to get the
1998 # resulting matrix. [3], verified against [4] (page 308).
1999 normalize = 1 / nsamples_treatments + 1 / nsamples_treatments[None].T
2001 # the standard error is used in the computation of the tukey criterion and
2002 # finding the p-values.
2003 stand_err = np.sqrt(normalize * mse / 2)
2005 # the mean difference is the test statistic.
2006 mean_differences = means[None].T - means
2008 # Calculate the t-statistic to use within the survival function of the
2009 # studentized range to get the p-value.
2010 t_stat = np.abs(mean_differences) / stand_err
2012 params = t_stat, ntreatments, nobs - ntreatments
2013 pvalues = distributions.studentized_range.sf(*params)
2015 return TukeyHSDResult(mean_differences, pvalues, ntreatments,
2016 nobs, stand_err)