Coverage for /usr/lib/python3/dist-packages/scipy/stats/_multicomp.py: 30%
116 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from __future__ import annotations
3import warnings
4from dataclasses import dataclass, field
5from typing import TYPE_CHECKING
7import numpy as np
9from scipy import stats
10from scipy.optimize import minimize_scalar
11from scipy.stats._common import ConfidenceInterval
12from scipy.stats._qmc import check_random_state
13from scipy.stats._stats_py import _var
15if TYPE_CHECKING:
16 import numpy.typing as npt
17 from scipy._lib._util import DecimalNumber, SeedType
18 from typing import Literal, Sequence # noqa: UP035
21__all__ = [
22 'dunnett'
23]
26@dataclass
27class DunnettResult:
28 """Result object returned by `scipy.stats.dunnett`.
30 Attributes
31 ----------
32 statistic : float ndarray
33 The computed statistic of the test for each comparison. The element
34 at index ``i`` is the statistic for the comparison between
35 groups ``i`` and the control.
36 pvalue : float ndarray
37 The computed p-value of the test for each comparison. The element
38 at index ``i`` is the p-value for the comparison between
39 group ``i`` and the control.
40 """
41 statistic: np.ndarray
42 pvalue: np.ndarray
43 _alternative: Literal['two-sided', 'less', 'greater'] = field(repr=False)
44 _rho: np.ndarray = field(repr=False)
45 _df: int = field(repr=False)
46 _std: float = field(repr=False)
47 _mean_samples: np.ndarray = field(repr=False)
48 _mean_control: np.ndarray = field(repr=False)
49 _n_samples: np.ndarray = field(repr=False)
50 _n_control: int = field(repr=False)
51 _rng: SeedType = field(repr=False)
52 _ci: ConfidenceInterval | None = field(default=None, repr=False)
53 _ci_cl: DecimalNumber | None = field(default=None, repr=False)
55 def __str__(self):
56 # Note: `__str__` prints the confidence intervals from the most
57 # recent call to `confidence_interval`. If it has not been called,
58 # it will be called with the default CL of .95.
59 if self._ci is None:
60 self.confidence_interval(confidence_level=.95)
61 s = (
62 "Dunnett's test"
63 f" ({self._ci_cl*100:.1f}% Confidence Interval)\n"
64 "Comparison Statistic p-value Lower CI Upper CI\n"
65 )
66 for i in range(self.pvalue.size):
67 s += (f" (Sample {i} - Control) {self.statistic[i]:>10.3f}"
68 f"{self.pvalue[i]:>10.3f}"
69 f"{self._ci.low[i]:>10.3f}"
70 f"{self._ci.high[i]:>10.3f}\n")
72 return s
74 def _allowance(
75 self, confidence_level: DecimalNumber = 0.95, tol: DecimalNumber = 1e-3
76 ) -> float:
77 """Allowance.
79 It is the quantity to add/subtract from the observed difference
80 between the means of observed groups and the mean of the control
81 group. The result gives confidence limits.
83 Parameters
84 ----------
85 confidence_level : float, optional
86 Confidence level for the computed confidence interval.
87 Default is .95.
88 tol : float, optional
89 A tolerance for numerical optimization: the allowance will produce
90 a confidence within ``10*tol*(1 - confidence_level)`` of the
91 specified level, or a warning will be emitted. Tight tolerances
92 may be impractical due to noisy evaluation of the objective.
93 Default is 1e-3.
95 Returns
96 -------
97 allowance : float
98 Allowance around the mean.
99 """
100 alpha = 1 - confidence_level
102 def pvalue_from_stat(statistic):
103 statistic = np.array(statistic)
104 sf = _pvalue_dunnett(
105 rho=self._rho, df=self._df,
106 statistic=statistic, alternative=self._alternative,
107 rng=self._rng
108 )
109 return abs(sf - alpha)/alpha
111 # Evaluation of `pvalue_from_stat` is noisy due to the use of RQMC to
112 # evaluate `multivariate_t.cdf`. `minimize_scalar` is not designed
113 # to tolerate a noisy objective function and may fail to find the
114 # minimum accurately. We mitigate this possibility with the validation
115 # step below, but implementation of a noise-tolerant root finder or
116 # minimizer would be a welcome enhancement. See gh-18150.
117 res = minimize_scalar(pvalue_from_stat, method='brent', tol=tol)
118 critical_value = res.x
120 # validation
121 # tol*10 because tol=1e-3 means we tolerate a 1% change at most
122 if res.success is False or res.fun >= tol*10:
123 warnings.warn(
124 "Computation of the confidence interval did not converge to "
125 "the desired level. The confidence level corresponding with "
126 f"the returned interval is approximately {alpha*(1+res.fun)}.",
127 stacklevel=3
128 )
130 # From [1] p. 1101 between (1) and (3)
131 allowance = critical_value*self._std*np.sqrt(
132 1/self._n_samples + 1/self._n_control
133 )
134 return abs(allowance)
136 def confidence_interval(
137 self, confidence_level: DecimalNumber = 0.95
138 ) -> ConfidenceInterval:
139 """Compute the confidence interval for the specified confidence level.
141 Parameters
142 ----------
143 confidence_level : float, optional
144 Confidence level for the computed confidence interval.
145 Default is .95.
147 Returns
148 -------
149 ci : ``ConfidenceInterval`` object
150 The object has attributes ``low`` and ``high`` that hold the
151 lower and upper bounds of the confidence intervals for each
152 comparison. The high and low values are accessible for each
153 comparison at index ``i`` for each group ``i``.
155 """
156 # check to see if the supplied confidence level matches that of the
157 # previously computed CI.
158 if (self._ci is not None) and (confidence_level == self._ci_cl):
159 return self._ci
161 if not (0 < confidence_level < 1):
162 raise ValueError("Confidence level must be between 0 and 1.")
164 allowance = self._allowance(confidence_level=confidence_level)
165 diff_means = self._mean_samples - self._mean_control
167 low = diff_means-allowance
168 high = diff_means+allowance
170 if self._alternative == 'greater':
171 high = [np.inf] * len(diff_means)
172 elif self._alternative == 'less':
173 low = [-np.inf] * len(diff_means)
175 self._ci_cl = confidence_level
176 self._ci = ConfidenceInterval(
177 low=low,
178 high=high
179 )
180 return self._ci
183def dunnett(
184 *samples: npt.ArrayLike, # noqa: D417
185 control: npt.ArrayLike,
186 alternative: Literal['two-sided', 'less', 'greater'] = "two-sided",
187 random_state: SeedType = None
188) -> DunnettResult:
189 """Dunnett's test: multiple comparisons of means against a control group.
191 This is an implementation of Dunnett's original, single-step test as
192 described in [1]_.
194 Parameters
195 ----------
196 sample1, sample2, ... : 1D array_like
197 The sample measurements for each experimental group.
198 control : 1D array_like
199 The sample measurements for the control group.
200 alternative : {'two-sided', 'less', 'greater'}, optional
201 Defines the alternative hypothesis.
203 The null hypothesis is that the means of the distributions underlying
204 the samples and control are equal. The following alternative
205 hypotheses are available (default is 'two-sided'):
207 * 'two-sided': the means of the distributions underlying the samples
208 and control are unequal.
209 * 'less': the means of the distributions underlying the samples
210 are less than the mean of the distribution underlying the control.
211 * 'greater': the means of the distributions underlying the
212 samples are greater than the mean of the distribution underlying
213 the control.
214 random_state : {None, int, `numpy.random.Generator`}, optional
215 If `random_state` is an int or None, a new `numpy.random.Generator` is
216 created using ``np.random.default_rng(random_state)``.
217 If `random_state` is already a ``Generator`` instance, then the
218 provided instance is used.
220 The random number generator is used to control the randomized
221 Quasi-Monte Carlo integration of the multivariate-t distribution.
223 Returns
224 -------
225 res : `~scipy.stats._result_classes.DunnettResult`
226 An object containing attributes:
228 statistic : float ndarray
229 The computed statistic of the test for each comparison. The element
230 at index ``i`` is the statistic for the comparison between
231 groups ``i`` and the control.
232 pvalue : float ndarray
233 The computed p-value of the test for each comparison. The element
234 at index ``i`` is the p-value for the comparison between
235 group ``i`` and the control.
237 And the following method:
239 confidence_interval(confidence_level=0.95) :
240 Compute the difference in means of the groups
241 with the control +- the allowance.
243 See Also
244 --------
245 tukey_hsd : performs pairwise comparison of means.
247 Notes
248 -----
249 Like the independent-sample t-test, Dunnett's test [1]_ is used to make
250 inferences about the means of distributions from which samples were drawn.
251 However, when multiple t-tests are performed at a fixed significance level,
252 the "family-wise error rate" - the probability of incorrectly rejecting the
253 null hypothesis in at least one test - will exceed the significance level.
254 Dunnett's test is designed to perform multiple comparisons while
255 controlling the family-wise error rate.
257 Dunnett's test compares the means of multiple experimental groups
258 against a single control group. Tukey's Honestly Significant Difference Test
259 is another multiple-comparison test that controls the family-wise error
260 rate, but `tukey_hsd` performs *all* pairwise comparisons between groups.
261 When pairwise comparisons between experimental groups are not needed,
262 Dunnett's test is preferable due to its higher power.
265 The use of this test relies on several assumptions.
267 1. The observations are independent within and among groups.
268 2. The observations within each group are normally distributed.
269 3. The distributions from which the samples are drawn have the same finite
270 variance.
272 References
273 ----------
274 .. [1] Charles W. Dunnett. "A Multiple Comparison Procedure for Comparing
275 Several Treatments with a Control."
276 Journal of the American Statistical Association, 50:272, 1096-1121,
277 :doi:`10.1080/01621459.1955.10501294`, 1955.
279 Examples
280 --------
281 In [1]_, the influence of drugs on blood count measurements on three groups
282 of animal is investigated.
284 The following table summarizes the results of the experiment in which
285 two groups received different drugs, and one group acted as a control.
286 Blood counts (in millions of cells per cubic millimeter) were recorded::
288 >>> import numpy as np
289 >>> control = np.array([7.40, 8.50, 7.20, 8.24, 9.84, 8.32])
290 >>> drug_a = np.array([9.76, 8.80, 7.68, 9.36])
291 >>> drug_b = np.array([12.80, 9.68, 12.16, 9.20, 10.55])
293 We would like to see if the means between any of the groups are
294 significantly different. First, visually examine a box and whisker plot.
296 >>> import matplotlib.pyplot as plt
297 >>> fig, ax = plt.subplots(1, 1)
298 >>> ax.boxplot([control, drug_a, drug_b])
299 >>> ax.set_xticklabels(["Control", "Drug A", "Drug B"]) # doctest: +SKIP
300 >>> ax.set_ylabel("mean") # doctest: +SKIP
301 >>> plt.show()
303 Note the overlapping interquartile ranges of the drug A group and control
304 group and the apparent separation between the drug B group and control
305 group.
307 Next, we will use Dunnett's test to assess whether the difference
308 between group means is significant while controlling the family-wise error
309 rate: the probability of making any false discoveries.
310 Let the null hypothesis be that the experimental groups have the same
311 mean as the control and the alternative be that an experimental group does
312 not have the same mean as the control. We will consider a 5% family-wise
313 error rate to be acceptable, and therefore we choose 0.05 as the threshold
314 for significance.
316 >>> from scipy.stats import dunnett
317 >>> res = dunnett(drug_a, drug_b, control=control)
318 >>> res.pvalue
319 array([0.62004941, 0.0059035 ]) # may vary
321 The p-value corresponding with the comparison between group A and control
322 exceeds 0.05, so we do not reject the null hypothesis for that comparison.
323 However, the p-value corresponding with the comparison between group B
324 and control is less than 0.05, so we consider the experimental results
325 to be evidence against the null hypothesis in favor of the alternative:
326 group B has a different mean than the control group.
328 """
329 samples_, control_, rng = _iv_dunnett(
330 samples=samples, control=control,
331 alternative=alternative, random_state=random_state
332 )
334 rho, df, n_group, n_samples, n_control = _params_dunnett(
335 samples=samples_, control=control_
336 )
338 statistic, std, mean_control, mean_samples = _statistic_dunnett(
339 samples_, control_, df, n_samples, n_control
340 )
342 pvalue = _pvalue_dunnett(
343 rho=rho, df=df, statistic=statistic, alternative=alternative, rng=rng
344 )
346 return DunnettResult(
347 statistic=statistic, pvalue=pvalue,
348 _alternative=alternative,
349 _rho=rho, _df=df, _std=std,
350 _mean_samples=mean_samples,
351 _mean_control=mean_control,
352 _n_samples=n_samples,
353 _n_control=n_control,
354 _rng=rng
355 )
358def _iv_dunnett(
359 samples: Sequence[npt.ArrayLike],
360 control: npt.ArrayLike,
361 alternative: Literal['two-sided', 'less', 'greater'],
362 random_state: SeedType
363) -> tuple[list[np.ndarray], np.ndarray, SeedType]:
364 """Input validation for Dunnett's test."""
365 rng = check_random_state(random_state)
367 if alternative not in {'two-sided', 'less', 'greater'}:
368 raise ValueError(
369 "alternative must be 'less', 'greater' or 'two-sided'"
370 )
372 ndim_msg = "Control and samples groups must be 1D arrays"
373 n_obs_msg = "Control and samples groups must have at least 1 observation"
375 control = np.asarray(control)
376 samples_ = [np.asarray(sample) for sample in samples]
378 # samples checks
379 samples_control: list[np.ndarray] = samples_ + [control]
380 for sample in samples_control:
381 if sample.ndim > 1:
382 raise ValueError(ndim_msg)
384 if sample.size < 1:
385 raise ValueError(n_obs_msg)
387 return samples_, control, rng
390def _params_dunnett(
391 samples: list[np.ndarray], control: np.ndarray
392) -> tuple[np.ndarray, int, int, np.ndarray, int]:
393 """Specific parameters for Dunnett's test.
395 Degree of freedom is the number of observations minus the number of groups
396 including the control.
397 """
398 n_samples = np.array([sample.size for sample in samples])
400 # From [1] p. 1100 d.f. = (sum N)-(p+1)
401 n_sample = n_samples.sum()
402 n_control = control.size
403 n = n_sample + n_control
404 n_groups = len(samples)
405 df = n - n_groups - 1
407 # From [1] p. 1103 rho_ij = 1/sqrt((N0/Ni+1)(N0/Nj+1))
408 rho = n_control/n_samples + 1
409 rho = 1/np.sqrt(rho[:, None] * rho[None, :])
410 np.fill_diagonal(rho, 1)
412 return rho, df, n_groups, n_samples, n_control
415def _statistic_dunnett(
416 samples: list[np.ndarray], control: np.ndarray, df: int,
417 n_samples: np.ndarray, n_control: int
418) -> tuple[np.ndarray, float, np.ndarray, np.ndarray]:
419 """Statistic of Dunnett's test.
421 Computation based on the original single-step test from [1].
422 """
423 mean_control = np.mean(control)
424 mean_samples = np.array([np.mean(sample) for sample in samples])
425 all_samples = [control] + samples
426 all_means = np.concatenate([[mean_control], mean_samples])
428 # Variance estimate s^2 from [1] Eq. 1
429 s2 = np.sum([_var(sample, mean=mean)*sample.size
430 for sample, mean in zip(all_samples, all_means)]) / df
431 std = np.sqrt(s2)
433 # z score inferred from [1] unlabeled equation after Eq. 1
434 z = (mean_samples - mean_control) / np.sqrt(1/n_samples + 1/n_control)
436 return z / std, std, mean_control, mean_samples
439def _pvalue_dunnett(
440 rho: np.ndarray, df: int, statistic: np.ndarray,
441 alternative: Literal['two-sided', 'less', 'greater'],
442 rng: SeedType = None
443) -> np.ndarray:
444 """pvalue from the multivariate t-distribution.
446 Critical values come from the multivariate student-t distribution.
447 """
448 statistic = statistic.reshape(-1, 1)
450 mvt = stats.multivariate_t(shape=rho, df=df, seed=rng)
451 if alternative == "two-sided":
452 statistic = abs(statistic)
453 pvalue = 1 - mvt.cdf(statistic, lower_limit=-statistic)
454 elif alternative == "greater":
455 pvalue = 1 - mvt.cdf(statistic, lower_limit=-np.inf)
456 else:
457 pvalue = 1 - mvt.cdf(np.inf, lower_limit=statistic)
459 return np.atleast_1d(pvalue)