Byar's Confidence Interval Calculation

Question

I was looking for advice on confidence intervals. I am currently working on automating the calculation of some of our indicators using Byar's method (outlined in Formula's 4a and 4b on page 7 of this document). I have coded these formulas out in python and I'm getting to-be-expected results with the following functions:

z = 1.96
def byars_lower(count, denominator, rate):
    c = 1 / (9 * count)
    b = 3 * sqrt(count)
    lower_o = count * ((1 - c - (z / b)) ** 3)
    lower_ci = (lower_o / denominator) * rate
    return lower_ci

def byars_higher(count, denominator, rate):
    c = 1 / (9 * (count + 1))
    b = 3 * (sqrt(count) + 1)
    upper_o = (count + 1) * ((1 - c + (z / b)) ** 3)
    upper_ci = (upper_o / denominator) * rate
    return upper_ci

The current team practice is to use an excel add-in to generate the CIs that basically have the following formulas:

lower CI = =IF(A2=0,0,IF(A2<389,CHIINV(0.5+95/200,2*A2)/2,A2*(1-1/(9*A2)-NORMSINV(0.5+95/200)/3/SQRT(A2))^3))/F2*C2

upper CI = =IF(A2<389,CHIINV(0.5-95/200,2*A2+2)/2,(A2+1)*(1-1/(9*(A2+1))+NORMSINV(0.5+95/200)/3/SQRT(A2+1))^3)/F2*C2

where: A2 = numerator, F2 = denominator & C2 = rate

There is a small discrepancy between these two methods that I can't reconcile (0.0005% - 0.001%), and they need to match to pass quality assurance. Would it be due to the CHIINV function in the excel formulas that would cause the difference?

Answer 1

So, after some more in depth analysis, there were a few issues:

z = 1.96 is stated in the documentation, but its actual value should have been z = ndtri(1 - alpha / 2) using from scipy.special import ndtri.

Line 3 in the byars_higher function should read b = 3 * (sqrt(count + 1)).

This solves all of the discrepancies for numerators above 389. I had to put a conditional statement in to deal with the numbers under 389. The chi2.ppf function from scipy.stats import chi2 replicates CHIINV from excel. The if statements are:

def byars_lower(count, denominator, rate):
    if count < 389:
        b = (chi2.ppf((alpha * 2), (count * 2)) / 2)
        c = b / denominator
        lower_ci = c * rate
        return lower_ci

def byars_higher(count, denominator, rate):
    if count < 389:
        b = chi2.ppf(1 - (alpha / 2), 2 * count + 2) / 2
        c = b / denominator
        upper_ci = c * rate
        return upper_ci

It now replicates the excel function above.