Coverage for /usr/lib/python3/dist-packages/scipy/stats/_relative_risk.py: 26%
57 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import operator
2from dataclasses import dataclass
3import numpy as np
4from scipy.special import ndtri
5from ._common import ConfidenceInterval
8def _validate_int(n, bound, name):
9 msg = f'{name} must be an integer not less than {bound}, but got {n!r}'
10 try:
11 n = operator.index(n)
12 except TypeError:
13 raise TypeError(msg) from None
14 if n < bound:
15 raise ValueError(msg)
16 return n
19@dataclass
20class RelativeRiskResult:
21 """
22 Result of `scipy.stats.contingency.relative_risk`.
24 Attributes
25 ----------
26 relative_risk : float
27 This is::
29 (exposed_cases/exposed_total) / (control_cases/control_total)
31 exposed_cases : int
32 The number of "cases" (i.e. occurrence of disease or other event
33 of interest) among the sample of "exposed" individuals.
34 exposed_total : int
35 The total number of "exposed" individuals in the sample.
36 control_cases : int
37 The number of "cases" among the sample of "control" or non-exposed
38 individuals.
39 control_total : int
40 The total number of "control" individuals in the sample.
42 Methods
43 -------
44 confidence_interval :
45 Compute the confidence interval for the relative risk estimate.
46 """
48 relative_risk: float
49 exposed_cases: int
50 exposed_total: int
51 control_cases: int
52 control_total: int
54 def confidence_interval(self, confidence_level=0.95):
55 """
56 Compute the confidence interval for the relative risk.
58 The confidence interval is computed using the Katz method
59 (i.e. "Method C" of [1]_; see also [2]_, section 3.1.2).
61 Parameters
62 ----------
63 confidence_level : float, optional
64 The confidence level to use for the confidence interval.
65 Default is 0.95.
67 Returns
68 -------
69 ci : ConfidenceInterval instance
70 The return value is an object with attributes ``low`` and
71 ``high`` that hold the confidence interval.
73 References
74 ----------
75 .. [1] D. Katz, J. Baptista, S. P. Azen and M. C. Pike, "Obtaining
76 confidence intervals for the risk ratio in cohort studies",
77 Biometrics, 34, 469-474 (1978).
78 .. [2] Hardeo Sahai and Anwer Khurshid, Statistics in Epidemiology,
79 CRC Press LLC, Boca Raton, FL, USA (1996).
82 Examples
83 --------
84 >>> from scipy.stats.contingency import relative_risk
85 >>> result = relative_risk(exposed_cases=10, exposed_total=75,
86 ... control_cases=12, control_total=225)
87 >>> result.relative_risk
88 2.5
89 >>> result.confidence_interval()
90 ConfidenceInterval(low=1.1261564003469628, high=5.549850800541033)
91 """
92 if not 0 <= confidence_level <= 1:
93 raise ValueError('confidence_level must be in the interval '
94 '[0, 1].')
96 # Handle edge cases where either exposed_cases or control_cases
97 # is zero. We follow the convention of the R function riskratio
98 # from the epitools library.
99 if self.exposed_cases == 0 and self.control_cases == 0:
100 # relative risk is nan.
101 return ConfidenceInterval(low=np.nan, high=np.nan)
102 elif self.exposed_cases == 0:
103 # relative risk is 0.
104 return ConfidenceInterval(low=0.0, high=np.nan)
105 elif self.control_cases == 0:
106 # relative risk is inf
107 return ConfidenceInterval(low=np.nan, high=np.inf)
109 alpha = 1 - confidence_level
110 z = ndtri(1 - alpha/2)
111 rr = self.relative_risk
113 # Estimate of the variance of log(rr) is
114 # var(log(rr)) = 1/exposed_cases - 1/exposed_total +
115 # 1/control_cases - 1/control_total
116 # and the standard error is the square root of that.
117 se = np.sqrt(1/self.exposed_cases - 1/self.exposed_total +
118 1/self.control_cases - 1/self.control_total)
119 delta = z*se
120 katz_lo = rr*np.exp(-delta)
121 katz_hi = rr*np.exp(delta)
122 return ConfidenceInterval(low=katz_lo, high=katz_hi)
125def relative_risk(exposed_cases, exposed_total, control_cases, control_total):
126 """
127 Compute the relative risk (also known as the risk ratio).
129 This function computes the relative risk associated with a 2x2
130 contingency table ([1]_, section 2.2.3; [2]_, section 3.1.2). Instead
131 of accepting a table as an argument, the individual numbers that are
132 used to compute the relative risk are given as separate parameters.
133 This is to avoid the ambiguity of which row or column of the contingency
134 table corresponds to the "exposed" cases and which corresponds to the
135 "control" cases. Unlike, say, the odds ratio, the relative risk is not
136 invariant under an interchange of the rows or columns.
138 Parameters
139 ----------
140 exposed_cases : nonnegative int
141 The number of "cases" (i.e. occurrence of disease or other event
142 of interest) among the sample of "exposed" individuals.
143 exposed_total : positive int
144 The total number of "exposed" individuals in the sample.
145 control_cases : nonnegative int
146 The number of "cases" among the sample of "control" or non-exposed
147 individuals.
148 control_total : positive int
149 The total number of "control" individuals in the sample.
151 Returns
152 -------
153 result : instance of `~scipy.stats._result_classes.RelativeRiskResult`
154 The object has the float attribute ``relative_risk``, which is::
156 rr = (exposed_cases/exposed_total) / (control_cases/control_total)
158 The object also has the method ``confidence_interval`` to compute
159 the confidence interval of the relative risk for a given confidence
160 level.
162 See Also
163 --------
164 odds_ratio
166 Notes
167 -----
168 The R package epitools has the function `riskratio`, which accepts
169 a table with the following layout::
171 disease=0 disease=1
172 exposed=0 (ref) n00 n01
173 exposed=1 n10 n11
175 With a 2x2 table in the above format, the estimate of the CI is
176 computed by `riskratio` when the argument method="wald" is given,
177 or with the function `riskratio.wald`.
179 For example, in a test of the incidence of lung cancer among a
180 sample of smokers and nonsmokers, the "exposed" category would
181 correspond to "is a smoker" and the "disease" category would
182 correspond to "has or had lung cancer".
184 To pass the same data to ``relative_risk``, use::
186 relative_risk(n11, n10 + n11, n01, n00 + n01)
188 .. versionadded:: 1.7.0
190 References
191 ----------
192 .. [1] Alan Agresti, An Introduction to Categorical Data Analysis
193 (second edition), Wiley, Hoboken, NJ, USA (2007).
194 .. [2] Hardeo Sahai and Anwer Khurshid, Statistics in Epidemiology,
195 CRC Press LLC, Boca Raton, FL, USA (1996).
197 Examples
198 --------
199 >>> from scipy.stats.contingency import relative_risk
201 This example is from Example 3.1 of [2]_. The results of a heart
202 disease study are summarized in the following table::
204 High CAT Low CAT Total
205 -------- ------- -----
206 CHD 27 44 71
207 No CHD 95 443 538
209 Total 122 487 609
211 CHD is coronary heart disease, and CAT refers to the level of
212 circulating catecholamine. CAT is the "exposure" variable, and
213 high CAT is the "exposed" category. So the data from the table
214 to be passed to ``relative_risk`` is::
216 exposed_cases = 27
217 exposed_total = 122
218 control_cases = 44
219 control_total = 487
221 >>> result = relative_risk(27, 122, 44, 487)
222 >>> result.relative_risk
223 2.4495156482861398
225 Find the confidence interval for the relative risk.
227 >>> result.confidence_interval(confidence_level=0.95)
228 ConfidenceInterval(low=1.5836990926700116, high=3.7886786315466354)
230 The interval does not contain 1, so the data supports the statement
231 that high CAT is associated with greater risk of CHD.
232 """
233 # Relative risk is a trivial calculation. The nontrivial part is in the
234 # `confidence_interval` method of the RelativeRiskResult class.
236 exposed_cases = _validate_int(exposed_cases, 0, "exposed_cases")
237 exposed_total = _validate_int(exposed_total, 1, "exposed_total")
238 control_cases = _validate_int(control_cases, 0, "control_cases")
239 control_total = _validate_int(control_total, 1, "control_total")
241 if exposed_cases > exposed_total:
242 raise ValueError('exposed_cases must not exceed exposed_total.')
243 if control_cases > control_total:
244 raise ValueError('control_cases must not exceed control_total.')
246 if exposed_cases == 0 and control_cases == 0:
247 # relative risk is 0/0.
248 rr = np.nan
249 elif exposed_cases == 0:
250 # relative risk is 0/nonzero
251 rr = 0.0
252 elif control_cases == 0:
253 # relative risk is nonzero/0.
254 rr = np.inf
255 else:
256 p1 = exposed_cases / exposed_total
257 p2 = control_cases / control_total
258 rr = p1 / p2
259 return RelativeRiskResult(relative_risk=rr,
260 exposed_cases=exposed_cases,
261 exposed_total=exposed_total,
262 control_cases=control_cases,
263 control_total=control_total)