`statspai.epi`¶

epi ¶

Epidemiology domain primitives (sp.epi).

Fills the gap the article calls out — statspai already has the heavy epidemiological causal machinery (IPW, G-formula, MSM, target trial), but lacked the entry-level statistical primitives that clinicians, epidemiologists, and public-health researchers reach for first.

Modelled after R's epiR, epitools, and fmsb.

import statspai as sp sp.epi.odds_ratio(50, 20, 30, 40) sp.epi.relative_risk(50, 950, 10, 990) sp.epi.mantel_haenszel(tables_2x2xK) sp.epi.direct_standardize(events, pop, standard_weights) sp.epi.bradford_hill(strength=1.0, temporality=1.0, consistency=0.5, ...)

OR2x2Result `dataclass` ¶

Bases: ResultProtocolMixin

Result of a 2x2 odds-ratio calculation.

RR2x2Result `dataclass` ¶

Bases: ResultProtocolMixin

Result of a 2x2 relative-risk (risk-ratio) calculation.

RD2x2Result `dataclass` ¶

Bases: ResultProtocolMixin

Result of a 2x2 risk-difference calculation.

ARResult `dataclass` ¶

Bases: ResultProtocolMixin

Attributable-risk quantities (Levin 1953, Miettinen 1974).

IRRResult `dataclass` ¶

Bases: ResultProtocolMixin

Incidence rate ratio from person-time data.

DiagnosticTestResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for binary diagnostic-test performance metrics.

Returned by :func:sensitivity_specificity / :func:diagnostic_test. Holds sensitivity and specificity (with Wilson-score CIs), predictive values (ppv, npv), likelihood ratios (lr_pos, lr_neg), prevalence, and the raw confusion cells tp/fp/fn/tn.

Examples:

>>> import statspai as sp
>>> res = sp.sensitivity_specificity(tp=90, fn=10, fp=5, tn=95)
>>> isinstance(res, sp.DiagnosticTestResult)
True
>>> res.sensitivity
0.9
>>> res.tp
90

ROCResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for ROC-curve coordinates and AUC inference.

Returned by :func:roc_curve. Holds the sweep thresholds and the corresponding true/false positive rates (tpr, fpr), plus the auc with its Hanley-McNeil standard error (auc_se) and CI (auc_ci).

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> y = np.array([0, 0, 0, 1, 1, 1])
>>> s = np.array([0.1, 0.2, 0.3, 0.7, 0.8, 0.9])  # perfectly separable
>>> roc = sp.roc_curve(y, s)
>>> isinstance(roc, sp.ROCResult)
True
>>> roc.auc
1.0

KappaResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for Cohen's kappa inter-rater agreement.

Returned by :func:cohen_kappa. Holds kappa with its standard error and CI, the observed and expected agreement, the number of n_categories, the weights scheme, and the z/p inference. The :meth:interpretation method maps kappa to a Landis-Koch label.

Examples:

>>> import statspai as sp
>>> a = [0, 1, 2, 0, 1, 2]
>>> b = [0, 1, 2, 0, 1, 2]  # perfect agreement
>>> k = sp.cohen_kappa(a, b)
>>> isinstance(k, sp.KappaResult)
True
>>> k.kappa
1.0
>>> k.interpretation()
'almost perfect agreement'

odds_ratio ¶

odds_ratio(a: Any, b: Optional[float] = None, c: Optional[float] = None, d: Optional[float] = None, *, method: str = 'woolf', alpha: float = 0.05) -> OR2x2Result

Odds ratio from a 2x2 table.

The standard epidemiology 2x2 layout is::

                Outcome+   Outcome-
Exposed            a          b
Unexposed          c          d

Parameters:

Name	Type	Description	Default
`a`	`float`	Cell counts, or pass a 2x2 array-like as `a`.	required
`b`	`float`	Cell counts, or pass a 2x2 array-like as `a`.	required
`c`	`float`	Cell counts, or pass a 2x2 array-like as `a`.	required
`d`	`float`	Cell counts, or pass a 2x2 array-like as `a`.	required
`method`	`('woolf', 'exact')`	Confidence-interval method. "woolf" uses the asymptotic log-OR standard error; "exact" uses the Fisher-style conditional non-central hypergeometric CI (via :func:`scipy.stats.fisher_exact`).	`"woolf"`
`alpha`	`float`		`0.05`

Returns:

Type	Description
`OR2x2Result`

Examples:

>>> import statspai as sp
>>> res = sp.epi.odds_ratio(50, 20, 30, 40)
>>> round(res.estimate, 3)
3.333

relative_risk ¶

relative_risk(a: Any, b: Optional[float] = None, c: Optional[float] = None, d: Optional[float] = None, *, alpha: float = 0.05) -> RR2x2Result

Relative risk (risk ratio) with Katz log-RR confidence interval.

Uses the Haldane correction when any cell is zero.

Examples:

>>> import statspai as sp
>>> res = sp.epi.relative_risk(40, 60, 20, 80)
>>> round(res.estimate, 3)
2.0
>>> round(res.risk_exposed, 3), round(res.risk_unexposed, 3)
(0.4, 0.2)

risk_difference ¶

risk_difference(a: Any, b: Optional[float] = None, c: Optional[float] = None, d: Optional[float] = None, *, method: str = 'wald', alpha: float = 0.05) -> RD2x2Result

Risk difference with Wald or Newcombe CI.

Parameters:

Name	Type	Description	Default
`method`	`('wald', 'newcombe')`	Newcombe's hybrid score CI avoids the Wald overshoot problem near 0 or 1.	`"wald"`

Examples:

>>> import statspai as sp
>>> res = sp.epi.risk_difference(40, 60, 20, 80)
>>> round(res.estimate, 3)
0.2
>>> round(res.risk_exposed, 3), round(res.risk_unexposed, 3)
(0.4, 0.2)

attributable_risk ¶

attributable_risk(a: Any, b: Optional[float] = None, c: Optional[float] = None, d: Optional[float] = None, *, alpha: float = 0.05) -> ARResult

Attributable fractions in exposed + in population (Levin 1953).

Computes: - AF_exposed = (RR - 1) / RR - PAF = P_e * (RR - 1) / [1 + P_e * (RR - 1)]

where P_e is prevalence of exposure. CI for PAF uses the delta method on log(1 - PAF).

Examples:

>>> import statspai as sp
>>> res = sp.epi.attributable_risk(40, 60, 20, 80)
>>> round(res.ar_exposed, 3)
0.5
>>> round(res.paf, 4), round(res.prevalence_exposed, 3)
(0.3333, 0.5)

incidence_rate_ratio ¶

incidence_rate_ratio(events_exposed: float, pt_exposed: float, events_unexposed: float, pt_unexposed: float, *, alpha: float = 0.05, method: str = 'exact') -> IRRResult

Person-time incidence rate ratio with exact Poisson CI.

Parameters:

Name	Type	Description	Default
`events_exposed`	`float`	Event counts.	required
`events_unexposed`	`float`	Event counts.	required
`pt_exposed`	`float`	Person-time at risk in each group (any time unit, as long as consistent).	required
`pt_unexposed`	`float`	Person-time at risk in each group (any time unit, as long as consistent).	required
`method`	`('exact', 'wald')`	"exact" uses the F-distribution-based Poisson CI (Breslow-Day); "wald" uses log-rate SE.	`"exact"`

Examples:

>>> import statspai as sp
>>> res = sp.epi.incidence_rate_ratio(30, 1000.0, 15, 1000.0)
>>> round(res.estimate, 3)
2.0
>>> round(res.rate_exposed, 3), round(res.rate_unexposed, 3)
(0.03, 0.015)

number_needed_to_treat ¶

number_needed_to_treat(a: Any, b: Optional[float] = None, c: Optional[float] = None, d: Optional[float] = None, *, alpha: float = 0.05) -> NNTResult

Number needed to treat (or harm), defined as |1 / RD|.

Propagates the Wald CI for RD. Interpretation convention: negative RD -> NNT-Benefit (treatment reduces risk); positive RD -> NNT-Harm.

Examples:

>>> import statspai as sp
>>> res = sp.epi.number_needed_to_treat(20, 80, 40, 60)
>>> round(res.estimate, 2)
5.0
>>> round(res.risk_difference, 3)  # treated arm has lower risk
-0.2

prevalence_ratio ¶

prevalence_ratio(*args: Any, **kwargs: Any) -> RR2x2Result

Prevalence ratio (cross-sectional RR); mathematically identical to :func:relative_risk when called on a 2x2 prevalence table. Distinguished for semantic clarity in cross-sectional studies.

Examples:

>>> import statspai as sp
>>> res = sp.epi.prevalence_ratio(40, 60, 20, 80)
>>> round(res.estimate, 3)
2.0
>>> res.method
'prevalence-ratio'

References

zou2004modified

mantel_haenszel ¶

mantel_haenszel(tables: Union[Sequence, ndarray], *, measure: str = 'OR', alpha: float = 0.05) -> MantelHaenszelResult

Mantel-Haenszel pooled OR / RR across K strata.

Parameters:

Name	Type	Description	Default
`tables`	array-like, shape ``(K, 2, 2)``	Each stratum's 2x2 table with layout `[[a_k, b_k], [c_k, d_k]]` (exposure x outcome).	required
`measure`	`('OR', 'RR')`	Pooled measure. Use :func:`mantel_haenszel_rate` for person-time IRR.	`"OR"`
`alpha`	`float`	Two-sided CI level.	`0.05`

Returns:

Type	Description
`MantelHaenszelResult`

Examples:

Pool the exposure-outcome odds ratio across two confounder strata, each a 2x2 table [[exposed-case, exposed-noncase], [unexposed-case, unexposed-noncase]]:

>>> import statspai as sp
>>> tables = [
...     [[10, 90], [5, 95]],
...     [[20, 80], [12, 88]],
... ]
>>> res = sp.mantel_haenszel(tables, measure="OR")
>>> type(res).__name__
'MantelHaenszelResult'
>>> res.measure
'OR'
>>> res.n_strata
2
>>> bool(res.estimate > 1.0)  # exposure raises the odds in both strata
True
>>> bool(res.ci[0] <= res.estimate <= res.ci[1])
True

References

mantel1959statistical

breslow_day_test ¶

breslow_day_test(tables: Union[Sequence, ndarray], *, tarone_correction: bool = True) -> tuple[float, float]

Breslow-Day test for homogeneity of the odds ratio across strata.

Parameters:

Name	Type	Description	Default
`tables`	array-like, shape ``(K, 2, 2)``		required
`tarone_correction`	`bool`	Apply Tarone's correction (recommended; Tarone 1985).	`True`

Returns:

Name	Type	Description
`chi2`	`float`
`p_value`	`float`

Examples:

Test whether the odds ratio is constant across two strata. Here the per-stratum ORs are close, so the test does not reject homogeneity:

>>> import statspai as sp
>>> tables = [
...     [[10, 90], [5, 95]],
...     [[20, 80], [12, 88]],
... ]
>>> chi2, p = sp.breslow_day_test(tables)
>>> bool(chi2 >= 0.0)
True
>>> bool(0.0 <= p <= 1.0)
True
>>> bool(p > 0.05)  # fail to reject OR homogeneity across strata
True

References

breslow1980statistical tarone1985heterogeneity

direct_standardize ¶

direct_standardize(events: Sequence[float], population: Sequence[float], standard_weights: Sequence[float], *, alpha: float = 0.05) -> StandardizedRateResult

Direct standardization of a rate.

The standardized rate is:

r_std = sum_k (w_k * events_k / population_k)

where w_k are the relative weights of a standard population (they are normalized internally to sum to 1).

Parameters:

Name	Type	Description	Default
`events`	`array - like`	Event counts in each stratum of the study population.	required
`population`	`array - like`	Denominator (person-time or population size) in each stratum.	required
`standard_weights`	`array - like`	Standard population size or proportion per stratum. Will be normalized to sum to 1.	required
`alpha`	`float`		`0.05`

Returns:

Type	Description
`StandardizedRateResult`

Notes

SE is computed by the delta method on the weighted sum of stratum rates, treating events as Poisson.

Examples:

>>> import statspai as sp
>>> events = [30, 50, 80]          # event counts per age stratum
>>> population = [1000, 1000, 500]
>>> standard_weights = [4000, 4000, 2000]   # external standard population
>>> res = sp.direct_standardize(events, population, standard_weights)
>>> round(res.rate, 4)
0.064

indirect_standardize ¶

indirect_standardize(observed: float, events_reference: Sequence[float], population_reference: Sequence[float], population_study: Sequence[float], *, alpha: float = 0.05) -> SMRResult

Indirect standardization -> Standardized Morbidity/Mortality Ratio.

Expected events = sum_k (rate_ref_k * pop_study_k), where rate_ref_k = events_reference_k / population_reference_k. SMR = observed / expected.

CI uses exact Poisson (Byar's approximation / Garwood).

Examples:

>>> import statspai as sp
>>> res = sp.indirect_standardize(
...     observed=120,
...     events_reference=[30, 50, 80],
...     population_reference=[1000, 1000, 500],
...     population_study=[800, 1200, 600],
... )
>>> round(res.expected, 4)
180.0
>>> round(res.smr, 4)
0.6667

diagnostic_test ¶

diagnostic_test(*args: Any, **kwargs: Any) -> DiagnosticTestResult

Alias for :func:sensitivity_specificity.

Examples:

>>> import statspai as sp
>>> res = sp.diagnostic_test(tp=90, fn=10, fp=5, tn=95)
>>> res.sensitivity  # 90 / (90 + 10)
0.9
>>> res.specificity  # 95 / (95 + 5)
0.95

sensitivity_specificity ¶

sensitivity_specificity(y_true: Any = None, y_pred: Any = None, *, tp: Optional[int] = None, fn: Optional[int] = None, fp: Optional[int] = None, tn: Optional[int] = None, alpha: float = 0.05) -> DiagnosticTestResult

Sensitivity and specificity with Wilson-score CIs.

Parameters:

Name	Type	Description	Default
`y_true`	`array - like`	Reference and predicted binary labels (0/1).	`None`
`y_pred`	`array - like`	Reference and predicted binary labels (0/1).	`None`
`tp`	`int`	Pre-computed confusion cells. Use instead of `y_true`/ `y_pred` when you already have counts.	`None`
`fn`	`int`	Pre-computed confusion cells. Use instead of `y_true`/ `y_pred` when you already have counts.	`None`
`fp`	`int`	Pre-computed confusion cells. Use instead of `y_true`/ `y_pred` when you already have counts.	`None`
`tn`	`int`	Pre-computed confusion cells. Use instead of `y_true`/ `y_pred` when you already have counts.	`None`
`alpha`	`float`		`0.05`

Examples:

>>> import statspai as sp
>>> res = sp.sensitivity_specificity(tp=90, fn=10, fp=5, tn=95)
>>> res.sensitivity  # 90 / (90 + 10)
0.9
>>> res.specificity  # 95 / (95 + 5)
0.95
>>> res.ppv  # 90 / (90 + 5)
0.9473684210526315

References

[@altman1994statistics]

roc_curve ¶

roc_curve(y_true: Any, scores: Any, *, alpha: float = 0.05) -> ROCResult

ROC curve with Hanley-McNeil (1982) AUC standard error.

Parameters:

Name	Type	Description	Default
`y_true`	`array-like of {0, 1}`		required
`scores`	`array-like of continuous predictions (higher = more "positive")`		required

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> y = np.array([0, 0, 0, 1, 1, 1])
>>> s = np.array([0.1, 0.2, 0.3, 0.7, 0.8, 0.9])  # perfectly separable
>>> roc = sp.roc_curve(y, s)
>>> roc.auc
1.0
>>> bool(roc.auc_ci[0] <= roc.auc <= roc.auc_ci[1])
True

References

[@hanley1982meaning]

auc ¶

auc(y_true: Any, scores: Any) -> float

Shortcut: just return the AUC.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> y = np.array([0, 0, 0, 1, 1, 1])
>>> s = np.array([0.1, 0.2, 0.3, 0.7, 0.8, 0.9])  # perfectly separable
>>> sp.auc(y, s)
1.0

cohen_kappa ¶

cohen_kappa(rater_a: Any, rater_b: Any, *, weights: str = 'unweighted', alpha: float = 0.05) -> KappaResult

Cohen's (1960) kappa for two raters.

Parameters:

Name	Type	Description	Default
`rater_a`	`array - like`	Same-length sequences of category labels from two raters.	required
`rater_b`	`array - like`	Same-length sequences of category labels from two raters.	required
`weights`	`('unweighted', 'linear', 'quadratic')`	Weighting scheme for disagreements across an ordered category scale. "unweighted" recovers the classic Cohen kappa.	`"unweighted"`

Examples:

>>> import statspai as sp
>>> a = [0, 1, 2, 0, 1, 2]
>>> b = [0, 1, 2, 0, 1, 2]  # perfect agreement
>>> k = sp.cohen_kappa(a, b)
>>> k.kappa
1.0
>>> k.n_categories
3
>>> c = [0, 1, 2, 0, 2, 1]  # two disagreements vs a
>>> kc = sp.cohen_kappa(a, c)
>>> bool(kc.kappa < 1.0)
True

References

[@cohen1960coefficient]

statspai.epi¶

epi ¶

OR2x2Result dataclass ¶

RR2x2Result dataclass ¶

RD2x2Result dataclass ¶

ARResult dataclass ¶

IRRResult dataclass ¶

DiagnosticTestResult dataclass ¶

ROCResult dataclass ¶

KappaResult dataclass ¶

odds_ratio ¶

relative_risk ¶

risk_difference ¶

attributable_risk ¶

incidence_rate_ratio ¶

number_needed_to_treat ¶

prevalence_ratio ¶

mantel_haenszel ¶

breslow_day_test ¶

direct_standardize ¶

indirect_standardize ¶

diagnostic_test ¶

sensitivity_specificity ¶

roc_curve ¶

auc ¶

cohen_kappa ¶

`statspai.epi`¶

OR2x2Result `dataclass` ¶

RR2x2Result `dataclass` ¶

RD2x2Result `dataclass` ¶

ARResult `dataclass` ¶

IRRResult `dataclass` ¶

DiagnosticTestResult `dataclass` ¶

ROCResult `dataclass` ¶

KappaResult `dataclass` ¶