`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 round(sp.epi.odds_ratio(50, 20, 30, 40).estimate, 3) 3.333 round(sp.epi.relative_risk(50, 950, 10, 990).estimate, 3) 5.0 tables_2x2xK = [[[10, 20], [5, 25]], [[8, 15], [6, 30]]] sp.epi.mantel_haenszel(tables_2x2xK).n_strata 2 events, pop = [10, 20, 30], [1000, 2000, 1500] standard_weights = [0.3, 0.4, 0.3] round(sp.epi.direct_standardize(events, pop, standard_weights).rate, 4) 0.013 sp.epi.bradford_hill(strength=1.0, temporality=1.0, consistency=0.5).total 2.5

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 or exact 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

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'

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 standard error (auc_se; Hanley-McNeil or DeLong, see se_method) 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

ARResult `dataclass` ¶

Bases: ResultProtocolMixin

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

IRRResult `dataclass` ¶

Bases: ResultProtocolMixin

Incidence rate ratio from person-time data.

NNTResult `dataclass` ¶

Bases: ResultProtocolMixin

Number needed to treat (benefit) or harm.

estimate is 1 / |RD| and ci the interval for that magnitude. ci_signed holds the reciprocals 1 / RD of the risk-difference limits in ascending order (what R epiR reports). When the risk-difference interval contains zero, ci_disjoint is True: the compatible NNTs are 1 / |L| or more on the benefit side and 1 / U or more on the harm side (Altman's "NNTB x to infinity to NNTH y"), ci_signed holds (1/L, 1/U) and ci is (1 / max(|L|, U), inf).

OR2x2Result `dataclass` ¶

Bases: ResultProtocolMixin

Result of a 2x2 odds-ratio calculation.

RD2x2Result `dataclass` ¶

Bases: ResultProtocolMixin

Result of a 2x2 risk-difference calculation.

RR2x2Result `dataclass` ¶

Bases: ResultProtocolMixin

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

auc ¶

auc(y_true: Any, scores: Any, weights: Any = None) -> float

Shortcut: just return the AUC (weighted when weights is given).

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]

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

ndcg ¶

ndcg(y_true: Any, scores: Any, k: float = 0.01) -> float

Normalised discounted cumulative gain of the top-ranked cases.

For a rare outcome (accounting fraud, default) the question is not how well a score ranks everyone, which is what the AUC answers, but how many true cases sit at the very top of the list an examiner would work through. NDCG at k looks at the k highest scores only, credits a true case at rank i with 1 / log2(i + 1), and divides by the credit of a perfect ranking, so 1 means the top of the list is all true cases (or all of them, when there are fewer than k).

Parameters:

Name	Type	Description	Default
`y_true`	`array - like`	0/1 outcome.	required
`scores`	`array - like`	Predicted score; higher means more likely.	required
`k`	`float`	A fraction in (0, 1): the top `round(n * k)` cases, 0.01 being the top 1%. A number of 1 or more is the count itself and must be a whole number.	`0.01`

Returns:

Type	Description
`float`	NDCG at `k` in [0, 1]; 0 when there is no true case at all.

Raises:

Type	Description
`MethodIncompatibility`	(a `ValueError`.) `y_true` holds a value other than 0 and 1, the inputs differ in length or hold missing values, or `k` selects no case.

Notes

Tied scores share their credit: every rank a group of tied cases occupies earns the group's share of true cases, so the value is the average of NDCG over all orderings of the ties and does not depend on the order of the rows (sklearn.metrics.ndcg_score with its default ignore_ties=False). Without ties this is the plain ranked sum, and equals farr::ndcg, which breaks ties in data order. Until 1.40.2 ties were broken in data order here too.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> y = np.array([1, 0, 1, 0, 0, 0, 0, 0, 0, 0])
>>> s = np.array([.9, .8, .7, .6, .5, .4, .3, .2, .1, .0])
>>> round(sp.ndcg(y, s, k=3), 4)        # hits at ranks 1 and 3 of top 3
0.9197
>>> sp.ndcg(y, s, k=1)
1.0

References

jarvelin2002cumulated, bao2020detecting

roc_curve ¶

roc_curve(y_true: Any, scores: Any, *, alpha: float = 0.05, se_method: str = 'delong', weights: Any = None) -> ROCResult

ROC curve with the AUC and its 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
`alpha`	`float`		`0.05`
`weights`	`array - like`	Non-negative observation weights. The rates become weighted shares and the AUC the weighted Mann-Whitney probability `sum_ij w_i w_j [1(s_i > s_j) + 1(s_i = s_j)/2] / (W+ W-)`. The main use is as a balance check: with `y_true` the exposure, `scores` the propensity score and `weights` the propensity weights, an AUC near 0.5 says the score no longer separates the two groups in the weighted sample. No standard error is reported with weights (`auc_se` and `auc_ci` are NaN): the weights are themselves estimated from the scores, and the DeLong and Hanley variances do not account for that.	`None`
`se_method`	`('delong', 'hanley', 'hanley-empirical')`	`"hanley"`: the Hanley-McNeil (1982) variance with the exponential approximations `Q1 = A/(2-A)`, `Q2 = 2A^2/(1+A)`. `"hanley-empirical"`: the same variance with `Q1` and `Q2` estimated from the data (ties weighted 1/3), which is what Stata `roctab, hanley` reports. `"delong"`: the DeLong-DeLong- Clarke-Pearson placement-value variance (R `pROC::var` / `ci.auc` default and Stata `roctab`'s default). .. versionchanged:: 1.30.0 The default moved from `"hanley"` (an approximation no reference uses by default) to `"delong"`, the default of both `pROC` and `roctab`. Pass `se_method="hanley"` for the old standard error.	`"delong"`

Notes

The curve is traced over the distinct score values: thresholds holds them in decreasing order and tpr / fpr the rates of the rule "positive when score >= threshold". The AUC is therefore the Mann-Whitney probability P(S+ > S-) + P(S+ = S-)/2 exactly, with tied scores counted as half. The confidence interval is AUC +/- z se truncated to [0, 1] (as pROC::ci.auc does).

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]

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, ci_method: str = 'wilson') -> DiagnosticTestResult

Sensitivity and specificity with Wilson-score (or exact) 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`
`ci_method`	`('wilson', 'exact')`	Interval for sensitivity and specificity: Wilson score (R `epiR::epi.tests(method="wilson")`) or Clopper-Pearson exact (`epi.tests` default, Stata `diagt`).	`"wilson"`

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]

attributable_risk ¶

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

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

Computes, from the raw counts: - AF_exposed = (RR - 1) / RR - PAF = P_e * (RR - 1) / [1 + P_e * (RR - 1)] = (R - R_0) / R

where P_e is the proportion exposed, R the overall risk and R_0 the risk among the unexposed (Stata cs / csi).

The PAF interval is the delta method on log(1 - PAF), where 1 - PAF = R_0 / R. prevalence says whether the exposure prevalence P_e = (a + b) / N is an estimate or a design constant:

"estimated" (default): the table is one random sample (a cross-sectional survey, or a cohort enrolled without regard to exposure), so P_e is estimated from it and its sampling error enters the interval. The variance is the multinomial delta method on the four cells. This is Stata punaf, vce(unconditional) after glm y x, family(binomial) link(log) vce(robust) except for the robust N / (N - 1) factor, which is left out here.
"fixed": the exposure groups were sampled separately with sizes fixed by design and P_e is known. The variance is the Katz variance of log RR scaled by d log(1 - PAF) / d log RR, as Stata punaf with the default conditional margins variance.

Under random sampling the "fixed" interval is too narrow when the outcome is common: with P_e = 0.3, R_0 = 0.45, RR = 2 and N = 400 it covered 91.4% in a simulation against 94.7% for "estimated".

The interval needs events in both exposure groups. With a = 0 or c = 0 the point estimates are still the plug-in values (PAF = 1 when no unexposed case occurred) but paf_ci is (nan, nan) with a warning; Stata cs reports no interval there either.

.. versionchanged:: 1.42 prevalence="estimated" is the default; the interval used to hold P_e fixed in every case ("fixed" keeps it).

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": the exact conditional (Clopper-Pearson) interval and the exact conditional binomial test, as R `poisson.test(c(D1, D0), c(PT1, PT0))`; with non-integer counts the p-value is the Wald one. "wald": `log IRR +/- z se` with `se = sqrt(1/D1 + 1/D0)` and the Wald test, a zero count replaced by 0.5 in both.	`"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.

When the RD interval (L, U) excludes zero, ci is (1/|U|, 1/|L|) sorted, matching R epiR::epi.2by2. When it includes zero the set of compatible NNTs is disjoint: 1/|L| or more on one side and 1/U or more on the other, reported by Altman (1998) as "NNTB 1/|L| to infinity to NNTH 1/U". ci_disjoint is then True, ci_signed is (1/L, 1/U) (epiR's pair) and ci is the interval for the magnitude, (1 / max(|L|, U), inf).

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

References

Altman, D. G. (1998). Confidence intervals for the number needed to treat. BMJ 317(7168), 1309-1312. doi:10.1136/bmj.317.7168.1309

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, correction: Optional[str] = None) -> 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')`	Inference method. "woolf" uses the asymptotic log-OR standard error and the Wald test; "exact" uses the exact conditional (noncentral hypergeometric) interval and Fisher's exact test. The point estimate `estimate` is the sample OR `ad / (bc)` under both methods, as Stata `cc` / `cci`.	`"woolf"`
`alpha`	`float`		`0.05`
`correction`	`(None, 'haldane')`	Only used by `method="woolf"`. `None`: with a zero cell the estimate is the plain ratio (0 or `inf`; NaN for 0/0) and the Woolf interval and Wald p-value are NaN, with a warning, as Stata `cci, woolf` reports no interval there. `"haldane"`: add 0.5 to every cell when any cell is zero and report that ratio and its finite interval. .. versionchanged:: 1.42 The Haldane correction used to be applied by default, so a zero cell moved the point estimate (`0, 10, 5, 20` gave 0.1774 instead of 0).	`None`

Returns:

Type	Description
`OR2x2Result`	With `method="exact"`, `conditional_mle` also holds the conditional maximum-likelihood OR, the estimate R `fisher.test` reports and the one the exact interval is built around.

Notes

The exact limits solve the two tail equations of the noncentral hypergeometric distribution to machine precision. R fisher.test and Stata cci, exact solve the same equations with looser root finders, so their limits differ from these in the fourth or fifth significant digit (R up to 1.2e-3 relative, Stata up to 2.2e-4 on the tables in the tests); R's conditional MLE differs by about 1e-5.

Examples:

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

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

relative_risk ¶

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

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

Parameters:

Name	Type	Description	Default
`a`	`float`	Cell counts (exposed cases, exposed non-cases, unexposed cases, unexposed non-cases), or pass a 2x2 array-like as `a`.	required
`b`	`float`	Cell counts (exposed cases, exposed non-cases, unexposed cases, unexposed non-cases), or pass a 2x2 array-like as `a`.	required
`c`	`float`	Cell counts (exposed cases, exposed non-cases, unexposed cases, unexposed non-cases), or pass a 2x2 array-like as `a`.	required
`d`	`float`	Cell counts (exposed cases, exposed non-cases, unexposed cases, unexposed non-cases), or pass a 2x2 array-like as `a`.	required
`alpha`	`float`		`0.05`
`correction`	`(None, 'haldane')`	`None`: the plain ratio `(a / (a + b)) / (c / (c + d))`, as Stata `cs` / `csi`. A zero `a` gives 0 and a zero `c` gives `inf`; the Katz interval and the Wald p-value need `a > 0` and `c > 0` and are NaN otherwise, with a warning (Stata reports no interval there either). `"haldane"`: add 0.5 to every cell when any cell is zero (Haldane-Anscombe) and report that ratio and its finite interval. .. versionchanged:: 1.42 The Haldane correction used to be applied by default, so a zero cell moved the point estimate (`0, 80, 20, 80` gave 0.0301 instead of 0).	`None`

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)

direct_standardize ¶

direct_standardize(events: Sequence[float], population: Sequence[float], standard_weights: Sequence[float], *, alpha: float = 0.05, ci_method: str = 'gamma', variance: str = 'poisson') -> 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`
`ci_method`	`('gamma', 'lognormal', 'normal')`	`"lognormal"`: `exp(log r +/- z se / r)`. `"gamma"`: the Fay-Feuer gamma interval, as R `epitools::ageadjust.direct`. `"normal"`: `r +/- z se` with the lower bound truncated at 0, as Stata `dstdize` (use with `variance="binomial"`). .. versionchanged:: 1.30.0 The default moved from `"lognormal"` (no package reference) to `"gamma"`, the default of `epitools::ageadjust.direct`; with the default `variance="poisson"` the rate and interval equal epitools. Pass `ci_method="lognormal"` for the old interval.	`"gamma"`
`variance`	`('poisson', 'binomial')`	`"poisson"`: `sum w_k^2 events_k / population_k^2` (epitools). `"binomial"`: `sum w_k^2 r_k (1 - r_k) / population_k` (Stata `dstdize`).	`"poisson"`

Returns:

Type	Description
`StandardizedRateResult`

Notes

The standard error is the square root of the chosen variance of the weighted sum of stratum rates; .se stores it.

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, ci_method: str = 'exact') -> 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_method="exact" (default) is the exact Poisson (Garwood) interval for the observed count divided by the expected count, as Stata istdize; "lognormal" is SMR exp(+/- z / sqrt(O)), as R epitools::ageadjust.indirect. p_value is the exact two-sided Poisson test of SMR = 1: twice the smaller exact tail, capped at 1.

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

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

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

statspai.epi¶

epi ¶

DiagnosticTestResult dataclass ¶

KappaResult dataclass ¶

ROCResult dataclass ¶

ARResult dataclass ¶

IRRResult dataclass ¶

NNTResult dataclass ¶

OR2x2Result dataclass ¶

RD2x2Result dataclass ¶

RR2x2Result dataclass ¶

auc ¶

cohen_kappa ¶

diagnostic_test ¶

ndcg ¶

roc_curve ¶

sensitivity_specificity ¶

attributable_risk ¶

incidence_rate_ratio ¶

number_needed_to_treat ¶

odds_ratio ¶

prevalence_ratio ¶

relative_risk ¶

risk_difference ¶

direct_standardize ¶

indirect_standardize ¶

breslow_day_test ¶

mantel_haenszel ¶

`statspai.epi`¶

DiagnosticTestResult `dataclass` ¶

KappaResult `dataclass` ¶

ROCResult `dataclass` ¶

ARResult `dataclass` ¶

IRRResult `dataclass` ¶

NNTResult `dataclass` ¶

OR2x2Result `dataclass` ¶

RD2x2Result `dataclass` ¶

RR2x2Result `dataclass` ¶