Coverage for /usr/lib/python3/dist-packages/scipy/stats/_survival.py: 22%
173 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from __future__ import annotations
3from dataclasses import dataclass, field
4from typing import TYPE_CHECKING
5import warnings
7import numpy as np
8from scipy import special, interpolate, stats
9from scipy.stats._censored_data import CensoredData
10from scipy.stats._common import ConfidenceInterval
12if TYPE_CHECKING:
13 from typing import Literal
14 import numpy.typing as npt
17__all__ = ['ecdf', 'logrank']
20@dataclass
21class EmpiricalDistributionFunction:
22 """An empirical distribution function produced by `scipy.stats.ecdf`
24 Attributes
25 ----------
26 quantiles : ndarray
27 The unique values of the sample from which the
28 `EmpiricalDistributionFunction` was estimated.
29 probabilities : ndarray
30 The point estimates of the cumulative distribution function (CDF) or
31 its complement, the survival function (SF), corresponding with
32 `quantiles`.
33 """
34 quantiles: np.ndarray
35 probabilities: np.ndarray
36 # Exclude these from __str__
37 _n: np.ndarray = field(repr=False) # number "at risk"
38 _d: np.ndarray = field(repr=False) # number of "deaths"
39 _sf: np.ndarray = field(repr=False) # survival function for var estimate
40 _kind: str = field(repr=False) # type of function: "cdf" or "sf"
42 def __init__(self, q, p, n, d, kind):
43 self.probabilities = p
44 self.quantiles = q
45 self._n = n
46 self._d = d
47 self._sf = p if kind == 'sf' else 1 - p
48 self._kind = kind
50 f0 = 1 if kind == 'sf' else 0 # leftmost function value
51 f1 = 1 - f0
52 # fill_value can't handle edge cases at infinity
53 x = np.insert(q, [0, len(q)], [-np.inf, np.inf])
54 y = np.insert(p, [0, len(p)], [f0, f1])
55 # `or` conditions handle the case of empty x, points
56 self._f = interpolate.interp1d(x, y, kind='previous',
57 assume_sorted=True)
59 def evaluate(self, x):
60 """Evaluate the empirical CDF/SF function at the input.
62 Parameters
63 ----------
64 x : ndarray
65 Argument to the CDF/SF
67 Returns
68 -------
69 y : ndarray
70 The CDF/SF evaluated at the input
71 """
72 return self._f(x)
74 def plot(self, ax=None, **matplotlib_kwargs):
75 """Plot the empirical distribution function
77 Available only if ``matplotlib`` is installed.
79 Parameters
80 ----------
81 ax : matplotlib.axes.Axes
82 Axes object to draw the plot onto, otherwise uses the current Axes.
84 **matplotlib_kwargs : dict, optional
85 Keyword arguments passed directly to `matplotlib.axes.Axes.step`.
86 Unless overridden, ``where='post'``.
88 Returns
89 -------
90 lines : list of `matplotlib.lines.Line2D`
91 Objects representing the plotted data
92 """
93 try:
94 import matplotlib # noqa
95 except ModuleNotFoundError as exc:
96 message = "matplotlib must be installed to use method `plot`."
97 raise ModuleNotFoundError(message) from exc
99 if ax is None:
100 import matplotlib.pyplot as plt
101 ax = plt.gca()
103 kwargs = {'where': 'post'}
104 kwargs.update(matplotlib_kwargs)
106 delta = np.ptp(self.quantiles)*0.05 # how far past sample edge to plot
107 q = self.quantiles
108 q = [q[0] - delta] + list(q) + [q[-1] + delta]
110 return ax.step(q, self.evaluate(q), **kwargs)
112 def confidence_interval(self, confidence_level=0.95, *, method='linear'):
113 """Compute a confidence interval around the CDF/SF point estimate
115 Parameters
116 ----------
117 confidence_level : float, default: 0.95
118 Confidence level for the computed confidence interval
120 method : str, {"linear", "log-log"}
121 Method used to compute the confidence interval. Options are
122 "linear" for the conventional Greenwood confidence interval
123 (default) and "log-log" for the "exponential Greenwood",
124 log-negative-log-transformed confidence interval.
126 Returns
127 -------
128 ci : ``ConfidenceInterval``
129 An object with attributes ``low`` and ``high``, instances of
130 `~scipy.stats._result_classes.EmpiricalDistributionFunction` that
131 represent the lower and upper bounds (respectively) of the
132 confidence interval.
134 Notes
135 -----
136 Confidence intervals are computed according to the Greenwood formula
137 (``method='linear'``) or the more recent "exponential Greenwood"
138 formula (``method='log-log'``) as described in [1]_. The conventional
139 Greenwood formula can result in lower confidence limits less than 0
140 and upper confidence limits greater than 1; these are clipped to the
141 unit interval. NaNs may be produced by either method; these are
142 features of the formulas.
144 References
145 ----------
146 .. [1] Sawyer, Stanley. "The Greenwood and Exponential Greenwood
147 Confidence Intervals in Survival Analysis."
148 https://www.math.wustl.edu/~sawyer/handouts/greenwood.pdf
150 """
151 message = ("Confidence interval bounds do not implement a "
152 "`confidence_interval` method.")
153 if self._n is None:
154 raise NotImplementedError(message)
156 methods = {'linear': self._linear_ci,
157 'log-log': self._loglog_ci}
159 message = f"`method` must be one of {set(methods)}."
160 if method.lower() not in methods:
161 raise ValueError(message)
163 message = "`confidence_level` must be a scalar between 0 and 1."
164 confidence_level = np.asarray(confidence_level)[()]
165 if confidence_level.shape or not (0 <= confidence_level <= 1):
166 raise ValueError(message)
168 method_fun = methods[method.lower()]
169 low, high = method_fun(confidence_level)
171 message = ("The confidence interval is undefined at some observations."
172 " This is a feature of the mathematical formula used, not"
173 " an error in its implementation.")
174 if np.any(np.isnan(low) | np.isnan(high)):
175 warnings.warn(message, RuntimeWarning, stacklevel=2)
177 low, high = np.clip(low, 0, 1), np.clip(high, 0, 1)
178 low = EmpiricalDistributionFunction(self.quantiles, low, None, None,
179 self._kind)
180 high = EmpiricalDistributionFunction(self.quantiles, high, None, None,
181 self._kind)
182 return ConfidenceInterval(low, high)
184 def _linear_ci(self, confidence_level):
185 sf, d, n = self._sf, self._d, self._n
186 # When n == d, Greenwood's formula divides by zero.
187 # When s != 0, this can be ignored: var == inf, and CI is [0, 1]
188 # When s == 0, this results in NaNs. Produce an informative warning.
189 with np.errstate(divide='ignore', invalid='ignore'):
190 var = sf ** 2 * np.cumsum(d / (n * (n - d)))
192 se = np.sqrt(var)
193 z = special.ndtri(1 / 2 + confidence_level / 2)
195 z_se = z * se
196 low = self.probabilities - z_se
197 high = self.probabilities + z_se
199 return low, high
201 def _loglog_ci(self, confidence_level):
202 sf, d, n = self._sf, self._d, self._n
204 with np.errstate(divide='ignore', invalid='ignore'):
205 var = 1 / np.log(sf) ** 2 * np.cumsum(d / (n * (n - d)))
207 se = np.sqrt(var)
208 z = special.ndtri(1 / 2 + confidence_level / 2)
210 with np.errstate(divide='ignore'):
211 lnl_points = np.log(-np.log(sf))
213 z_se = z * se
214 low = np.exp(-np.exp(lnl_points + z_se))
215 high = np.exp(-np.exp(lnl_points - z_se))
216 if self._kind == "cdf":
217 low, high = 1-high, 1-low
219 return low, high
222@dataclass
223class ECDFResult:
224 """ Result object returned by `scipy.stats.ecdf`
226 Attributes
227 ----------
228 cdf : `~scipy.stats._result_classes.EmpiricalDistributionFunction`
229 An object representing the empirical cumulative distribution function.
230 sf : `~scipy.stats._result_classes.EmpiricalDistributionFunction`
231 An object representing the complement of the empirical cumulative
232 distribution function.
233 """
234 cdf: EmpiricalDistributionFunction
235 sf: EmpiricalDistributionFunction
237 def __init__(self, q, cdf, sf, n, d):
238 self.cdf = EmpiricalDistributionFunction(q, cdf, n, d, "cdf")
239 self.sf = EmpiricalDistributionFunction(q, sf, n, d, "sf")
242def _iv_CensoredData(
243 sample: npt.ArrayLike | CensoredData, param_name: str = 'sample'
244) -> CensoredData:
245 """Attempt to convert `sample` to `CensoredData`."""
246 if not isinstance(sample, CensoredData):
247 try: # takes care of input standardization/validation
248 sample = CensoredData(uncensored=sample)
249 except ValueError as e:
250 message = str(e).replace('uncensored', param_name)
251 raise type(e)(message) from e
252 return sample
255def ecdf(sample: npt.ArrayLike | CensoredData) -> ECDFResult:
256 """Empirical cumulative distribution function of a sample.
258 The empirical cumulative distribution function (ECDF) is a step function
259 estimate of the CDF of the distribution underlying a sample. This function
260 returns objects representing both the empirical distribution function and
261 its complement, the empirical survival function.
263 Parameters
264 ----------
265 sample : 1D array_like or `scipy.stats.CensoredData`
266 Besides array_like, instances of `scipy.stats.CensoredData` containing
267 uncensored and right-censored observations are supported. Currently,
268 other instances of `scipy.stats.CensoredData` will result in a
269 ``NotImplementedError``.
271 Returns
272 -------
273 res : `~scipy.stats._result_classes.ECDFResult`
274 An object with the following attributes.
276 cdf : `~scipy.stats._result_classes.EmpiricalDistributionFunction`
277 An object representing the empirical cumulative distribution
278 function.
279 sf : `~scipy.stats._result_classes.EmpiricalDistributionFunction`
280 An object representing the empirical survival function.
282 The `cdf` and `sf` attributes themselves have the following attributes.
284 quantiles : ndarray
285 The unique values in the sample that defines the empirical CDF/SF.
286 probabilities : ndarray
287 The point estimates of the probabilities corresponding with
288 `quantiles`.
290 And the following methods:
292 evaluate(x) :
293 Evaluate the CDF/SF at the argument.
295 plot(ax) :
296 Plot the CDF/SF on the provided axes.
298 confidence_interval(confidence_level=0.95) :
299 Compute the confidence interval around the CDF/SF at the values in
300 `quantiles`.
302 Notes
303 -----
304 When each observation of the sample is a precise measurement, the ECDF
305 steps up by ``1/len(sample)`` at each of the observations [1]_.
307 When observations are lower bounds, upper bounds, or both upper and lower
308 bounds, the data is said to be "censored", and `sample` may be provided as
309 an instance of `scipy.stats.CensoredData`.
311 For right-censored data, the ECDF is given by the Kaplan-Meier estimator
312 [2]_; other forms of censoring are not supported at this time.
314 Confidence intervals are computed according to the Greenwood formula or the
315 more recent "Exponential Greenwood" formula as described in [4]_.
317 References
318 ----------
319 .. [1] Conover, William Jay. Practical nonparametric statistics. Vol. 350.
320 John Wiley & Sons, 1999.
322 .. [2] Kaplan, Edward L., and Paul Meier. "Nonparametric estimation from
323 incomplete observations." Journal of the American statistical
324 association 53.282 (1958): 457-481.
326 .. [3] Goel, Manish Kumar, Pardeep Khanna, and Jugal Kishore.
327 "Understanding survival analysis: Kaplan-Meier estimate."
328 International journal of Ayurveda research 1.4 (2010): 274.
330 .. [4] Sawyer, Stanley. "The Greenwood and Exponential Greenwood Confidence
331 Intervals in Survival Analysis."
332 https://www.math.wustl.edu/~sawyer/handouts/greenwood.pdf
334 Examples
335 --------
336 **Uncensored Data**
338 As in the example from [1]_ page 79, five boys were selected at random from
339 those in a single high school. Their one-mile run times were recorded as
340 follows.
342 >>> sample = [6.23, 5.58, 7.06, 6.42, 5.20] # one-mile run times (minutes)
344 The empirical distribution function, which approximates the distribution
345 function of one-mile run times of the population from which the boys were
346 sampled, is calculated as follows.
348 >>> from scipy import stats
349 >>> res = stats.ecdf(sample)
350 >>> res.cdf.quantiles
351 array([5.2 , 5.58, 6.23, 6.42, 7.06])
352 >>> res.cdf.probabilities
353 array([0.2, 0.4, 0.6, 0.8, 1. ])
355 To plot the result as a step function:
357 >>> import matplotlib.pyplot as plt
358 >>> ax = plt.subplot()
359 >>> res.cdf.plot(ax)
360 >>> ax.set_xlabel('One-Mile Run Time (minutes)')
361 >>> ax.set_ylabel('Empirical CDF')
362 >>> plt.show()
364 **Right-censored Data**
366 As in the example from [1]_ page 91, the lives of ten car fanbelts were
367 tested. Five tests concluded because the fanbelt being tested broke, but
368 the remaining tests concluded for other reasons (e.g. the study ran out of
369 funding, but the fanbelt was still functional). The mileage driven
370 with the fanbelts were recorded as follows.
372 >>> broken = [77, 47, 81, 56, 80] # in thousands of miles driven
373 >>> unbroken = [62, 60, 43, 71, 37]
375 Precise survival times of the fanbelts that were still functional at the
376 end of the tests are unknown, but they are known to exceed the values
377 recorded in ``unbroken``. Therefore, these observations are said to be
378 "right-censored", and the data is represented using
379 `scipy.stats.CensoredData`.
381 >>> sample = stats.CensoredData(uncensored=broken, right=unbroken)
383 The empirical survival function is calculated as follows.
385 >>> res = stats.ecdf(sample)
386 >>> res.sf.quantiles
387 array([37., 43., 47., 56., 60., 62., 71., 77., 80., 81.])
388 >>> res.sf.probabilities
389 array([1. , 1. , 0.875, 0.75 , 0.75 , 0.75 , 0.75 , 0.5 , 0.25 , 0. ])
391 To plot the result as a step function:
393 >>> ax = plt.subplot()
394 >>> res.cdf.plot(ax)
395 >>> ax.set_xlabel('Fanbelt Survival Time (thousands of miles)')
396 >>> ax.set_ylabel('Empirical SF')
397 >>> plt.show()
399 """
400 sample = _iv_CensoredData(sample)
402 if sample.num_censored() == 0:
403 res = _ecdf_uncensored(sample._uncensor())
404 elif sample.num_censored() == sample._right.size:
405 res = _ecdf_right_censored(sample)
406 else:
407 # Support additional censoring options in follow-up PRs
408 message = ("Currently, only uncensored and right-censored data is "
409 "supported.")
410 raise NotImplementedError(message)
412 t, cdf, sf, n, d = res
413 return ECDFResult(t, cdf, sf, n, d)
416def _ecdf_uncensored(sample):
417 sample = np.sort(sample)
418 x, counts = np.unique(sample, return_counts=True)
420 # [1].81 "the fraction of [observations] that are less than or equal to x
421 events = np.cumsum(counts)
422 n = sample.size
423 cdf = events / n
425 # [1].89 "the relative frequency of the sample that exceeds x in value"
426 sf = 1 - cdf
428 at_risk = np.concatenate(([n], n - events[:-1]))
429 return x, cdf, sf, at_risk, counts
432def _ecdf_right_censored(sample):
433 # It is conventional to discuss right-censored data in terms of
434 # "survival time", "death", and "loss" (e.g. [2]). We'll use that
435 # terminology here.
436 # This implementation was influenced by the references cited and also
437 # https://www.youtube.com/watch?v=lxoWsVco_iM
438 # https://en.wikipedia.org/wiki/Kaplan%E2%80%93Meier_estimator
439 # In retrospect it is probably most easily compared against [3].
440 # Ultimately, the data needs to be sorted, so this implementation is
441 # written to avoid a separate call to `unique` after sorting. In hope of
442 # better performance on large datasets, it also computes survival
443 # probabilities at unique times only rather than at each observation.
444 tod = sample._uncensored # time of "death"
445 tol = sample._right # time of "loss"
446 times = np.concatenate((tod, tol))
447 died = np.asarray([1]*tod.size + [0]*tol.size)
449 # sort by times
450 i = np.argsort(times)
451 times = times[i]
452 died = died[i]
453 at_risk = np.arange(times.size, 0, -1)
455 # logical indices of unique times
456 j = np.diff(times, prepend=-np.inf, append=np.inf) > 0
457 j_l = j[:-1] # first instances of unique times
458 j_r = j[1:] # last instances of unique times
460 # get number at risk and deaths at each unique time
461 t = times[j_l] # unique times
462 n = at_risk[j_l] # number at risk at each unique time
463 cd = np.cumsum(died)[j_r] # cumulative deaths up to/including unique times
464 d = np.diff(cd, prepend=0) # deaths at each unique time
466 # compute survival function
467 sf = np.cumprod((n - d) / n)
468 cdf = 1 - sf
469 return t, cdf, sf, n, d
472@dataclass
473class LogRankResult:
474 """Result object returned by `scipy.stats.logrank`.
476 Attributes
477 ----------
478 statistic : float ndarray
479 The computed statistic (defined below). Its magnitude is the
480 square root of the magnitude returned by most other logrank test
481 implementations.
482 pvalue : float ndarray
483 The computed p-value of the test.
484 """
485 statistic: np.ndarray
486 pvalue: np.ndarray
489def logrank(
490 x: npt.ArrayLike | CensoredData,
491 y: npt.ArrayLike | CensoredData,
492 alternative: Literal['two-sided', 'less', 'greater'] = "two-sided"
493) -> LogRankResult:
494 r"""Compare the survival distributions of two samples via the logrank test.
496 Parameters
497 ----------
498 x, y : array_like or CensoredData
499 Samples to compare based on their empirical survival functions.
500 alternative : {'two-sided', 'less', 'greater'}, optional
501 Defines the alternative hypothesis.
503 The null hypothesis is that the survival distributions of the two
504 groups, say *X* and *Y*, are identical.
506 The following alternative hypotheses [4]_ are available (default is
507 'two-sided'):
509 * 'two-sided': the survival distributions of the two groups are not
510 identical.
511 * 'less': survival of group *X* is favored: the group *X* failure rate
512 function is less than the group *Y* failure rate function at some
513 times.
514 * 'greater': survival of group *Y* is favored: the group *X* failure
515 rate function is greater than the group *Y* failure rate function at
516 some times.
518 Returns
519 -------
520 res : `~scipy.stats._result_classes.LogRankResult`
521 An object containing attributes:
523 statistic : float ndarray
524 The computed statistic (defined below). Its magnitude is the
525 square root of the magnitude returned by most other logrank test
526 implementations.
527 pvalue : float ndarray
528 The computed p-value of the test.
530 See Also
531 --------
532 scipy.stats.ecdf
534 Notes
535 -----
536 The logrank test [1]_ compares the observed number of events to
537 the expected number of events under the null hypothesis that the two
538 samples were drawn from the same distribution. The statistic is
540 .. math::
542 Z_i = \frac{\sum_{j=1}^J(O_{i,j}-E_{i,j})}{\sqrt{\sum_{j=1}^J V_{i,j}}}
543 \rightarrow \mathcal{N}(0,1)
545 where
547 .. math::
549 E_{i,j} = O_j \frac{N_{i,j}}{N_j},
550 \qquad
551 V_{i,j} = E_{i,j} \left(\frac{N_j-O_j}{N_j}\right)
552 \left(\frac{N_j-N_{i,j}}{N_j-1}\right),
554 :math:`i` denotes the group (i.e. it may assume values :math:`x` or
555 :math:`y`, or it may be omitted to refer to the combined sample)
556 :math:`j` denotes the time (at which an event occured),
557 :math:`N` is the number of subjects at risk just before an event occured,
558 and :math:`O` is the observed number of events at that time.
560 The ``statistic`` :math:`Z_x` returned by `logrank` is the (signed) square
561 root of the statistic returned by many other implementations. Under the
562 null hypothesis, :math:`Z_x**2` is asymptotically distributed according to
563 the chi-squared distribution with one degree of freedom. Consequently,
564 :math:`Z_x` is asymptotically distributed according to the standard normal
565 distribution. The advantage of using :math:`Z_x` is that the sign
566 information (i.e. whether the observed number of events tends to be less
567 than or greater than the number expected under the null hypothesis) is
568 preserved, allowing `scipy.stats.logrank` to offer one-sided alternative
569 hypotheses.
571 References
572 ----------
573 .. [1] Mantel N. "Evaluation of survival data and two new rank order
574 statistics arising in its consideration."
575 Cancer Chemotherapy Reports, 50(3):163-170, PMID: 5910392, 1966
576 .. [2] Bland, Altman, "The logrank test", BMJ, 328:1073,
577 :doi:`10.1136/bmj.328.7447.1073`, 2004
578 .. [3] "Logrank test", Wikipedia,
579 https://en.wikipedia.org/wiki/Logrank_test
580 .. [4] Brown, Mark. "On the choice of variance for the log rank test."
581 Biometrika 71.1 (1984): 65-74.
582 .. [5] Klein, John P., and Melvin L. Moeschberger. Survival analysis:
583 techniques for censored and truncated data. Vol. 1230. New York:
584 Springer, 2003.
586 Examples
587 --------
588 Reference [2]_ compared the survival times of patients with two different
589 types of recurrent malignant gliomas. The samples below record the time
590 (number of weeks) for which each patient participated in the study. The
591 `scipy.stats.CensoredData` class is used because the data is
592 right-censored: the uncensored observations correspond with observed deaths
593 whereas the censored observations correspond with the patient leaving the
594 study for another reason.
596 >>> from scipy import stats
597 >>> x = stats.CensoredData(
598 ... uncensored=[6, 13, 21, 30, 37, 38, 49, 50,
599 ... 63, 79, 86, 98, 202, 219],
600 ... right=[31, 47, 80, 82, 82, 149]
601 ... )
602 >>> y = stats.CensoredData(
603 ... uncensored=[10, 10, 12, 13, 14, 15, 16, 17, 18, 20, 24, 24,
604 ... 25, 28,30, 33, 35, 37, 40, 40, 46, 48, 76, 81,
605 ... 82, 91, 112, 181],
606 ... right=[34, 40, 70]
607 ... )
609 We can calculate and visualize the empirical survival functions
610 of both groups as follows.
612 >>> import numpy as np
613 >>> import matplotlib.pyplot as plt
614 >>> ax = plt.subplot()
615 >>> ecdf_x = stats.ecdf(x)
616 >>> ecdf_x.sf.plot(ax, label='Astrocytoma')
617 >>> ecdf_y = stats.ecdf(y)
618 >>> ecdf_x.sf.plot(ax, label='Glioblastoma')
619 >>> ax.set_xlabel('Time to death (weeks)')
620 >>> ax.set_ylabel('Empirical SF')
621 >>> plt.legend()
622 >>> plt.show()
624 Visual inspection of the empirical survival functions suggests that the
625 survival times tend to be different between the two groups. To formally
626 assess whether the difference is significant at the 1% level, we use the
627 logrank test.
629 >>> res = stats.logrank(x=x, y=y)
630 >>> res.statistic
631 -2.73799...
632 >>> res.pvalue
633 0.00618...
635 The p-value is less than 1%, so we can consider the data to be evidence
636 against the null hypothesis in favor of the alternative that there is a
637 difference between the two survival functions.
639 """
640 # Input validation. `alternative` IV handled in `_normtest_finish` below.
641 x = _iv_CensoredData(sample=x, param_name='x')
642 y = _iv_CensoredData(sample=y, param_name='y')
644 # Combined sample. (Under H0, the two groups are identical.)
645 xy = CensoredData(
646 uncensored=np.concatenate((x._uncensored, y._uncensored)),
647 right=np.concatenate((x._right, y._right))
648 )
650 # Extract data from the combined sample
651 res = ecdf(xy)
652 idx = res.sf._d.astype(bool) # indices of observed events
653 times_xy = res.sf.quantiles[idx] # unique times of observed events
654 at_risk_xy = res.sf._n[idx] # combined number of subjects at risk
655 deaths_xy = res.sf._d[idx] # combined number of events
657 # Get the number at risk within each sample.
658 # First compute the number at risk in group X at each of the `times_xy`.
659 # Could use `interpolate_1d`, but this is more compact.
660 res_x = ecdf(x)
661 i = np.searchsorted(res_x.sf.quantiles, times_xy)
662 at_risk_x = np.append(res_x.sf._n, 0)[i] # 0 at risk after last time
663 # Subtract from the combined number at risk to get number at risk in Y
664 at_risk_y = at_risk_xy - at_risk_x
666 # Compute the variance.
667 num = at_risk_x * at_risk_y * deaths_xy * (at_risk_xy - deaths_xy)
668 den = at_risk_xy**2 * (at_risk_xy - 1)
669 # Note: when `at_risk_xy == 1`, we would have `at_risk_xy - 1 == 0` in the
670 # numerator and denominator. Simplifying the fraction symbolically, we
671 # would always find the overall quotient to be zero, so don't compute it.
672 i = at_risk_xy > 1
673 sum_var = np.sum(num[i]/den[i])
675 # Get the observed and expected number of deaths in group X
676 n_died_x = x._uncensored.size
677 sum_exp_deaths_x = np.sum(at_risk_x * (deaths_xy/at_risk_xy))
679 # Compute the statistic. This is the square root of that in references.
680 statistic = (n_died_x - sum_exp_deaths_x)/np.sqrt(sum_var)
682 # Equivalent to chi2(df=1).sf(statistic**2) when alternative='two-sided'
683 _, pvalue = stats._stats_py._normtest_finish(
684 z=statistic, alternative=alternative
685 )
687 return LogRankResult(statistic=statistic, pvalue=pvalue)