NUMERICAL AUDIT · 1 NOVEMBER 2026

statsmodels: the sqrtcc rate-ratio interval ignores exposures

confint_poisson_2indep(method="sqrtcc", compare="ratio") returns an interval for the ratio of counts rather than rates, so with unequal follow-up a nominal 95% interval covered the true ratio 46.6% of the time.

Xamit Kadirbekov, Daniyal Kadirbekov · AI-assisted investigation and writing

Checked 1 November 2026: PR #10660 is open, with no review yet; not merged or released. Local regression checks and limitations documented.

Upstream report #10660

Authors: Xamit Kadirbekov and Daniyal Kadirbekov (GERO). Published 1 November 2026.

Status: Checked 1 November 2026: PR #10660 is open, with no review yet; not merged or released. Local regression checks and limitations documented.

Published proposal

statsmodels PR #10660 — head commit 6be12a7c153021893567ea5679fee146893335cd, OPEN at submission on 1 November 2026.

What the function does

An event rate is the number of events divided by the time observed (the exposure), for example 41 infections in 2,969 patient-years. To compare two groups, trials and epidemiological studies use the rate ratio (x1/n1)/(x2/n2); dividing by exposure removes the effect of one group simply being observed longer. statsmodels.stats.rates.confint_poisson_2indep gives a confidence interval for this ratio by several methods. sqrtcc is the Price–Bonett interval based on square roots of the counts with a +0.5 correction. It is not the default (score is), but it is listed in the docstring as an available method.

The discrepancy

In stats/rates.py (lines 1870–1881 on main 5bb03e5; identical in release 0.15.0):

Exposures are ignored. The function returns an interval for the ratio of counts, not of rates. exposure1 and exposure2 are never used, while the other methods of the same function rescale by n2/n1.
The Fieller discriminant uses 0.25 * crit instead of 0.25 * crit**2, while the denominator, two lines below, uses crit**2. Its effect is small (about 0.3% in the bounds).

from statsmodels.stats.rates import confint_poisson_2indep as ci
ci(41, 2969, 15, 1106, method="score",  compare="ratio")   # (0.568, 1.824)   rate ratio = 1.018
ci(41, 2969, 15, 1106, method="mover",  compare="ratio")   # (0.564, 1.823)
ci(41, 2969, 15, 1106, method="sqrtcc", compare="ratio")   # (1.532, 5.067)   <- excludes 1
ci(41, 2969*1000, 15, 1106*1000, method="sqrtcc", compare="ratio")  # (1.532, 5.067), unchanged

The point estimate 1.018 lies outside the sqrtcc interval. The reference interval for this example is (0.573, 1.884). Setting both exposures to 1 also leaves the sqrtcc interval unchanged. With equal exposures (30/100 against 30/100) only the small second error remains: (0.5998, 1.6656) against the reference (0.6013, 1.6630).

Coverage of the nominal 95% interval in a simulation with true ratio 2, exposures 200 and 100, 2,000 replications: sqrtcc 46.6%, score 95.0%. In an earlier run with true ratio 1 and the same exposures (1,000 replications), sqrtcc covered 59.9% against 94.3% for the reference.

Reference and why this is not a convention

The reference is the Fieller interval for √(μ1/μ2) on the counts (Price & Bonett 2000, eq. 2.4, the equation the code comment cites), re-derived and evaluated in mpmath at 50 digits, then rescaled to rates by n2/n1. The other methods of the function agree with it. The function documents its target as the ratio of two Poisson rates and accepts exposures; an interval that does not move when one exposure is multiplied by 1,000 cannot be an interval for a rate ratio under any convention. We did not have access to the full text of the paper, so the check of the second error rests on the derivation, not on comparing the printed equation.

Fix and tests

Rescale the sqrtcc bounds by n2/n1 and use crit**2 in the discriminant.
Three new tests in stats/tests/test_rates_poisson.py: invariance to scaling both exposures, agreement with the reference on the example above, and sqrtcc against score on a balanced design. They fail before and pass after. The rest of the module: 1169 → 1172 passed. Three existing tests fail with NameError both before and after in the local environment.

Upstream status

As of 1 November 2026, PR #10660 is open, with no review or comment. Nothing has been merged or released.

Limits

The tests ran on an installed 0.15.0 wheel with the patched pure-Python module, since rates.py is identical on main 5bb03e5; this was not a fresh build of main. Users of the default score method are not affected. With x2 = 0 the upper bound should be infinite; it is meaningless both before and after the patch, and the PR does not address it. The method was added in #8166; we found no later issue or PR about sqrtcc, and #10354 and #10373 add input validation without touching this formula. This is a bounded GitHub search.

The example and the coverage simulation were re-run in a separate session on 1 November before submission. No reproduction archive is attached to this page yet.

Attribution and license

AI tools (Claude) helped find the cases, build the reference calculations and draft the patches and tests; the examples were re-run independently before submission. Xamit Kadirbekov and Daniyal Kadirbekov are the report authors and reviewed the submitted diff; upstream software remains credited to its developers. Report: CC BY 4.0. Upstream patches retain their upstream licenses.