0
$\begingroup$

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?

$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.