Coverage for /usr/lib/python3/dist-packages/scipy/stats/_odds_ratio.py: 15%
137 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import numpy as np
3from scipy.special import ndtri
4from scipy.optimize import brentq
5from ._discrete_distns import nchypergeom_fisher
6from ._common import ConfidenceInterval
9def _sample_odds_ratio(table):
10 """
11 Given a table [[a, b], [c, d]], compute a*d/(b*c).
13 Return nan if the numerator and denominator are 0.
14 Return inf if just the denominator is 0.
15 """
16 # table must be a 2x2 numpy array.
17 if table[1, 0] > 0 and table[0, 1] > 0:
18 oddsratio = table[0, 0] * table[1, 1] / (table[1, 0] * table[0, 1])
19 elif table[0, 0] == 0 or table[1, 1] == 0:
20 oddsratio = np.nan
21 else:
22 oddsratio = np.inf
23 return oddsratio
26def _solve(func):
27 """
28 Solve func(nc) = 0. func must be an increasing function.
29 """
30 # We could just as well call the variable `x` instead of `nc`, but we
31 # always call this function with functions for which nc (the noncentrality
32 # parameter) is the variable for which we are solving.
33 nc = 1.0
34 value = func(nc)
35 if value == 0:
36 return nc
38 # Multiplicative factor by which to increase or decrease nc when
39 # searching for a bracketing interval.
40 factor = 2.0
41 # Find a bracketing interval.
42 if value > 0:
43 nc /= factor
44 while func(nc) > 0:
45 nc /= factor
46 lo = nc
47 hi = factor*nc
48 else:
49 nc *= factor
50 while func(nc) < 0:
51 nc *= factor
52 lo = nc/factor
53 hi = nc
55 # lo and hi bracket the solution for nc.
56 nc = brentq(func, lo, hi, xtol=1e-13)
57 return nc
60def _nc_hypergeom_mean_inverse(x, M, n, N):
61 """
62 For the given noncentral hypergeometric parameters x, M, n,and N
63 (table[0,0], total, row 0 sum and column 0 sum, resp., of a 2x2
64 contingency table), find the noncentrality parameter of Fisher's
65 noncentral hypergeometric distribution whose mean is x.
66 """
67 nc = _solve(lambda nc: nchypergeom_fisher.mean(M, n, N, nc) - x)
68 return nc
71def _hypergeom_params_from_table(table):
72 # The notation M, n and N is consistent with stats.hypergeom and
73 # stats.nchypergeom_fisher.
74 x = table[0, 0]
75 M = table.sum()
76 n = table[0].sum()
77 N = table[:, 0].sum()
78 return x, M, n, N
81def _ci_upper(table, alpha):
82 """
83 Compute the upper end of the confidence interval.
84 """
85 if _sample_odds_ratio(table) == np.inf:
86 return np.inf
88 x, M, n, N = _hypergeom_params_from_table(table)
90 # nchypergeom_fisher.cdf is a decreasing function of nc, so we negate
91 # it in the lambda expression.
92 nc = _solve(lambda nc: -nchypergeom_fisher.cdf(x, M, n, N, nc) + alpha)
93 return nc
96def _ci_lower(table, alpha):
97 """
98 Compute the lower end of the confidence interval.
99 """
100 if _sample_odds_ratio(table) == 0:
101 return 0
103 x, M, n, N = _hypergeom_params_from_table(table)
105 nc = _solve(lambda nc: nchypergeom_fisher.sf(x - 1, M, n, N, nc) - alpha)
106 return nc
109def _conditional_oddsratio(table):
110 """
111 Conditional MLE of the odds ratio for the 2x2 contingency table.
112 """
113 x, M, n, N = _hypergeom_params_from_table(table)
114 # Get the bounds of the support. The support of the noncentral
115 # hypergeometric distribution with parameters M, n, and N is the same
116 # for all values of the noncentrality parameter, so we can use 1 here.
117 lo, hi = nchypergeom_fisher.support(M, n, N, 1)
119 # Check if x is at one of the extremes of the support. If so, we know
120 # the odds ratio is either 0 or inf.
121 if x == lo:
122 # x is at the low end of the support.
123 return 0
124 if x == hi:
125 # x is at the high end of the support.
126 return np.inf
128 nc = _nc_hypergeom_mean_inverse(x, M, n, N)
129 return nc
132def _conditional_oddsratio_ci(table, confidence_level=0.95,
133 alternative='two-sided'):
134 """
135 Conditional exact confidence interval for the odds ratio.
136 """
137 if alternative == 'two-sided':
138 alpha = 0.5*(1 - confidence_level)
139 lower = _ci_lower(table, alpha)
140 upper = _ci_upper(table, alpha)
141 elif alternative == 'less':
142 lower = 0.0
143 upper = _ci_upper(table, 1 - confidence_level)
144 else:
145 # alternative == 'greater'
146 lower = _ci_lower(table, 1 - confidence_level)
147 upper = np.inf
149 return lower, upper
152def _sample_odds_ratio_ci(table, confidence_level=0.95,
153 alternative='two-sided'):
154 oddsratio = _sample_odds_ratio(table)
155 log_or = np.log(oddsratio)
156 se = np.sqrt((1/table).sum())
157 if alternative == 'less':
158 z = ndtri(confidence_level)
159 loglow = -np.inf
160 loghigh = log_or + z*se
161 elif alternative == 'greater':
162 z = ndtri(confidence_level)
163 loglow = log_or - z*se
164 loghigh = np.inf
165 else:
166 # alternative is 'two-sided'
167 z = ndtri(0.5*confidence_level + 0.5)
168 loglow = log_or - z*se
169 loghigh = log_or + z*se
171 return np.exp(loglow), np.exp(loghigh)
174class OddsRatioResult:
175 """
176 Result of `scipy.stats.contingency.odds_ratio`. See the
177 docstring for `odds_ratio` for more details.
179 Attributes
180 ----------
181 statistic : float
182 The computed odds ratio.
184 * If `kind` is ``'sample'``, this is sample (or unconditional)
185 estimate, given by
186 ``table[0, 0]*table[1, 1]/(table[0, 1]*table[1, 0])``.
187 * If `kind` is ``'conditional'``, this is the conditional
188 maximum likelihood estimate for the odds ratio. It is
189 the noncentrality parameter of Fisher's noncentral
190 hypergeometric distribution with the same hypergeometric
191 parameters as `table` and whose mean is ``table[0, 0]``.
193 Methods
194 -------
195 confidence_interval :
196 Confidence interval for the odds ratio.
197 """
199 def __init__(self, _table, _kind, statistic):
200 # for now, no need to make _table and _kind public, since this sort of
201 # information is returned in very few `scipy.stats` results
202 self._table = _table
203 self._kind = _kind
204 self.statistic = statistic
206 def __repr__(self):
207 return f"OddsRatioResult(statistic={self.statistic})"
209 def confidence_interval(self, confidence_level=0.95,
210 alternative='two-sided'):
211 """
212 Confidence interval for the odds ratio.
214 Parameters
215 ----------
216 confidence_level: float
217 Desired confidence level for the confidence interval.
218 The value must be given as a fraction between 0 and 1.
219 Default is 0.95 (meaning 95%).
221 alternative : {'two-sided', 'less', 'greater'}, optional
222 The alternative hypothesis of the hypothesis test to which the
223 confidence interval corresponds. That is, suppose the null
224 hypothesis is that the true odds ratio equals ``OR`` and the
225 confidence interval is ``(low, high)``. Then the following options
226 for `alternative` are available (default is 'two-sided'):
228 * 'two-sided': the true odds ratio is not equal to ``OR``. There
229 is evidence against the null hypothesis at the chosen
230 `confidence_level` if ``high < OR`` or ``low > OR``.
231 * 'less': the true odds ratio is less than ``OR``. The ``low`` end
232 of the confidence interval is 0, and there is evidence against
233 the null hypothesis at the chosen `confidence_level` if
234 ``high < OR``.
235 * 'greater': the true odds ratio is greater than ``OR``. The
236 ``high`` end of the confidence interval is ``np.inf``, and there
237 is evidence against the null hypothesis at the chosen
238 `confidence_level` if ``low > OR``.
240 Returns
241 -------
242 ci : ``ConfidenceInterval`` instance
243 The confidence interval, represented as an object with
244 attributes ``low`` and ``high``.
246 Notes
247 -----
248 When `kind` is ``'conditional'``, the limits of the confidence
249 interval are the conditional "exact confidence limits" as described
250 by Fisher [1]_. The conditional odds ratio and confidence interval are
251 also discussed in Section 4.1.2 of the text by Sahai and Khurshid [2]_.
253 When `kind` is ``'sample'``, the confidence interval is computed
254 under the assumption that the logarithm of the odds ratio is normally
255 distributed with standard error given by::
257 se = sqrt(1/a + 1/b + 1/c + 1/d)
259 where ``a``, ``b``, ``c`` and ``d`` are the elements of the
260 contingency table. (See, for example, [2]_, section 3.1.3.2,
261 or [3]_, section 2.3.3).
263 References
264 ----------
265 .. [1] R. A. Fisher (1935), The logic of inductive inference,
266 Journal of the Royal Statistical Society, Vol. 98, No. 1,
267 pp. 39-82.
268 .. [2] H. Sahai and A. Khurshid (1996), Statistics in Epidemiology:
269 Methods, Techniques, and Applications, CRC Press LLC, Boca
270 Raton, Florida.
271 .. [3] Alan Agresti, An Introduction to Categorical Data Analyis
272 (second edition), Wiley, Hoboken, NJ, USA (2007).
273 """
274 if alternative not in ['two-sided', 'less', 'greater']:
275 raise ValueError("`alternative` must be 'two-sided', 'less' or "
276 "'greater'.")
278 if confidence_level < 0 or confidence_level > 1:
279 raise ValueError('confidence_level must be between 0 and 1')
281 if self._kind == 'conditional':
282 ci = self._conditional_odds_ratio_ci(confidence_level, alternative)
283 else:
284 ci = self._sample_odds_ratio_ci(confidence_level, alternative)
285 return ci
287 def _conditional_odds_ratio_ci(self, confidence_level=0.95,
288 alternative='two-sided'):
289 """
290 Confidence interval for the conditional odds ratio.
291 """
293 table = self._table
294 if 0 in table.sum(axis=0) or 0 in table.sum(axis=1):
295 # If both values in a row or column are zero, the p-value is 1,
296 # the odds ratio is NaN and the confidence interval is (0, inf).
297 ci = (0, np.inf)
298 else:
299 ci = _conditional_oddsratio_ci(table,
300 confidence_level=confidence_level,
301 alternative=alternative)
302 return ConfidenceInterval(low=ci[0], high=ci[1])
304 def _sample_odds_ratio_ci(self, confidence_level=0.95,
305 alternative='two-sided'):
306 """
307 Confidence interval for the sample odds ratio.
308 """
309 if confidence_level < 0 or confidence_level > 1:
310 raise ValueError('confidence_level must be between 0 and 1')
312 table = self._table
313 if 0 in table.sum(axis=0) or 0 in table.sum(axis=1):
314 # If both values in a row or column are zero, the p-value is 1,
315 # the odds ratio is NaN and the confidence interval is (0, inf).
316 ci = (0, np.inf)
317 else:
318 ci = _sample_odds_ratio_ci(table,
319 confidence_level=confidence_level,
320 alternative=alternative)
321 return ConfidenceInterval(low=ci[0], high=ci[1])
324def odds_ratio(table, *, kind='conditional'):
325 r"""
326 Compute the odds ratio for a 2x2 contingency table.
328 Parameters
329 ----------
330 table : array_like of ints
331 A 2x2 contingency table. Elements must be non-negative integers.
332 kind : str, optional
333 Which kind of odds ratio to compute, either the sample
334 odds ratio (``kind='sample'``) or the conditional odds ratio
335 (``kind='conditional'``). Default is ``'conditional'``.
337 Returns
338 -------
339 result : `~scipy.stats._result_classes.OddsRatioResult` instance
340 The returned object has two computed attributes:
342 statistic : float
343 * If `kind` is ``'sample'``, this is sample (or unconditional)
344 estimate, given by
345 ``table[0, 0]*table[1, 1]/(table[0, 1]*table[1, 0])``.
346 * If `kind` is ``'conditional'``, this is the conditional
347 maximum likelihood estimate for the odds ratio. It is
348 the noncentrality parameter of Fisher's noncentral
349 hypergeometric distribution with the same hypergeometric
350 parameters as `table` and whose mean is ``table[0, 0]``.
352 The object has the method `confidence_interval` that computes
353 the confidence interval of the odds ratio.
355 See Also
356 --------
357 scipy.stats.fisher_exact
358 relative_risk
360 Notes
361 -----
362 The conditional odds ratio was discussed by Fisher (see "Example 1"
363 of [1]_). Texts that cover the odds ratio include [2]_ and [3]_.
365 .. versionadded:: 1.10.0
367 References
368 ----------
369 .. [1] R. A. Fisher (1935), The logic of inductive inference,
370 Journal of the Royal Statistical Society, Vol. 98, No. 1,
371 pp. 39-82.
372 .. [2] Breslow NE, Day NE (1980). Statistical methods in cancer research.
373 Volume I - The analysis of case-control studies. IARC Sci Publ.
374 (32):5-338. PMID: 7216345. (See section 4.2.)
375 .. [3] H. Sahai and A. Khurshid (1996), Statistics in Epidemiology:
376 Methods, Techniques, and Applications, CRC Press LLC, Boca
377 Raton, Florida.
378 .. [4] Berger, Jeffrey S. et al. "Aspirin for the Primary Prevention of
379 Cardiovascular Events in Women and Men: A Sex-Specific
380 Meta-analysis of Randomized Controlled Trials."
381 JAMA, 295(3):306-313, :doi:`10.1001/jama.295.3.306`, 2006.
383 Examples
384 --------
385 In epidemiology, individuals are classified as "exposed" or
386 "unexposed" to some factor or treatment. If the occurrence of some
387 illness is under study, those who have the illness are often
388 classifed as "cases", and those without it are "noncases". The
389 counts of the occurrences of these classes gives a contingency
390 table::
392 exposed unexposed
393 cases a b
394 noncases c d
396 The sample odds ratio may be written ``(a/c) / (b/d)``. ``a/c`` can
397 be interpreted as the odds of a case occurring in the exposed group,
398 and ``b/d`` as the odds of a case occurring in the unexposed group.
399 The sample odds ratio is the ratio of these odds. If the odds ratio
400 is greater than 1, it suggests that there is a positive association
401 between being exposed and being a case.
403 Interchanging the rows or columns of the contingency table inverts
404 the odds ratio, so it is import to understand the meaning of labels
405 given to the rows and columns of the table when interpreting the
406 odds ratio.
408 In [4]_, the use of aspirin to prevent cardiovascular events in women
409 and men was investigated. The study notably concluded:
411 ...aspirin therapy reduced the risk of a composite of
412 cardiovascular events due to its effect on reducing the risk of
413 ischemic stroke in women [...]
415 The article lists studies of various cardiovascular events. Let's
416 focus on the ischemic stoke in women.
418 The following table summarizes the results of the experiment in which
419 participants took aspirin or a placebo on a regular basis for several
420 years. Cases of ischemic stroke were recorded::
422 Aspirin Control/Placebo
423 Ischemic stroke 176 230
424 No stroke 21035 21018
426 The question we ask is "Is there evidence that the aspirin reduces the
427 risk of ischemic stroke?"
429 Compute the odds ratio:
431 >>> from scipy.stats.contingency import odds_ratio
432 >>> res = odds_ratio([[176, 230], [21035, 21018]])
433 >>> res.statistic
434 0.7646037659999126
436 For this sample, the odds of getting an ischemic stroke for those who have
437 been taking aspirin are 0.76 times that of those
438 who have received the placebo.
440 To make statistical inferences about the population under study,
441 we can compute the 95% confidence interval for the odds ratio:
443 >>> res.confidence_interval(confidence_level=0.95)
444 ConfidenceInterval(low=0.6241234078749812, high=0.9354102892100372)
446 The 95% confidence interval for the conditional odds ratio is
447 approximately (0.62, 0.94).
449 The fact that the entire 95% confidence interval falls below 1 supports
450 the authors' conclusion that the aspirin was associated with a
451 statistically significant reduction in ischemic stroke.
452 """
453 if kind not in ['conditional', 'sample']:
454 raise ValueError("`kind` must be 'conditional' or 'sample'.")
456 c = np.asarray(table)
458 if c.shape != (2, 2):
459 raise ValueError(f"Invalid shape {c.shape}. The input `table` must be "
460 "of shape (2, 2).")
462 if not np.issubdtype(c.dtype, np.integer):
463 raise ValueError("`table` must be an array of integers, but got "
464 f"type {c.dtype}")
465 c = c.astype(np.int64)
467 if np.any(c < 0):
468 raise ValueError("All values in `table` must be nonnegative.")
470 if 0 in c.sum(axis=0) or 0 in c.sum(axis=1):
471 # If both values in a row or column are zero, the p-value is NaN and
472 # the odds ratio is NaN.
473 result = OddsRatioResult(_table=c, _kind=kind, statistic=np.nan)
474 return result
476 if kind == 'sample':
477 oddsratio = _sample_odds_ratio(c)
478 else: # kind is 'conditional'
479 oddsratio = _conditional_oddsratio(c)
481 result = OddsRatioResult(_table=c, _kind=kind, statistic=oddsratio)
482 return result