Formulas, Algorithms and Examples for Binomial Distributed Data Confidence Interval Calculation: Excess Risk, Relative Risk and Odds Ratio

Abstract

Medical studies often involve a comparison between two outcomes, each collected from a sample. The probability associated with, and confidence in the result of the study is of most importance, since one may argue that having been wrong with a percent could be what killed a patient. Sampling is usually done from a finite and discrete population and it follows a Bernoulli trial, leading to a contingency of two binomially distributed samples (better known as $2\times 2$ contingency table). Current guidelines recommend reporting relative measures of association (such as the relative risk and odds ratio) in conjunction with absolute measures of association (which include risk difference or excess risk). Because the distribution is discrete, the evaluation of the exact confidence interval for either of those measures of association is a mathematical challenge. Some alternate scenarios were analyzed (continuous vs. discrete; hypergeometric vs. binomial), and in the main case—bivariate binomial experiment—a strategy for providing exact p-values and confidence intervals is proposed. Algorithms implementing the strategy are given.

1. Introduction

Evaluation of the odds and the risks associated with diseases becomes more actual than never in the recent context of COVID-19 [1,2]. The proper evaluation of these parameters as well as their confidence intervals have a great impact [3]. The term odds denotes the probability that an event will occur divided by the probability that an event will not occur while the ‘Risk’ is the probability that an outcome will occur, expressed as the number of positive outcomes divided by the total number of outcomes; relative risk and odds ratios are simply the risk or odds, respectively, of an outcome in one group (exposed group, case group) divided by the risk or odds of the outcome in another group (control group, reference group). Current guidelines recommend reporting relative measures of association (such as the relative risk and odds ratio) in conjunction with absolute measures of association (which include risk difference or excess risk) [4]; Table 4 in [4] should be consulted for a detailed comparison between them.

There are no interpretations of these concepts that are at once simple, intuitive, correct, and foolproof; correct use and interpretation of these statistics requires an attention to detail which seems to tax the patience of working scientists [5]. A recent study noted that assessment of intervention effects in systematic reviews with meta-analysis deserves greater rigor with respect to insufficiently demonstrated traditional 95% confidence intervals and p-values [6] suggesting that the primary step of reporting the 95% confidence interval, and/or the p-value should be continued with other four steps (investigating statistical and clinical heterogeneity; dealing with problems of multiplicity due to multiple outcomes; trial sequential analysis; and Bayes factor analysis). In addition, to avoid erroneous interpretations, researchers should account for the potential impact of systematic errors (bias; on the meta-analysis results and publication bias due to statistical significance and clinical significance. The assessment of the clinical significance of an intervention effect should only be assessed if statistical significance has been obtained, but on the other hand, if statistical significance has been reached, then clinical significance must be assessed [6]. There is no doubt that in research practice and in training of researchers, we must engage appropriate methods including calculation of the confidence intervals within individual studies, and meta-analysis in the integration of multiple studies, in order to produce final conclusions [7].

Even if providing an accurately calculated confidence interval (CI) might strongly support the conclusions of a study, people are still reticent to use exact methods of calculating the CIs. Furthermore, some reputed statisticians suggest that “approximate is better than exact” [8]. To follow the same route, one could say “this concludes the needed proof”.

The aforementioned reticence is in the majority of cases due to the complexity of the calculation. For example, deciphering how Casella ([9] refining a previous work [10]) calculated the famous BSC interval for the binomial proportion is a challenge for someone not quite familiarized with mathematics, statistics and informatics. To briefly exemplify, one should assume that, generally, the people developing a method and the people meant to use it are not the same, and worse yet, may not have the same background.

New findings showed that the bias associated with the usual estimator of odds ratio and the relative risk may be substantial for rare events [11]. Further developments are expected from machine learning enabling the development of algorithms with minimal dependence or parametric assumptions [12].

2. Research Aim

When dealing with sampling, usually from a finite and discrete population, Value–Probability and Interval–Confidence facets for the $2\times 2$ contingency measures of association require a much better explanation in order to be thoroughly understood, and “at key solutions” need to be provided. Here, it is proposed that starting from the study design (Figure 1), first, to calculate the probability mass functions associated with the two binomial samples, second, to calculate the probability mass function of the case vs. control bivariate binomial experiment, and finally, to extract the parameters of interest (please note that most of the studies require calculation of at least two parameters). Cumulative distribution function and confidence intervals and/or p-values are calculated in a succession of two more steps (Figure 1); later, algorithms are given (Algorithms 1–5).

Algorithm 1: Probability mass function (PMF) for binomial distribution, Equation (9).

Input:z, o    //z←x or y, o←m or n

procedure PMF_B(z, o, &r)

If($z\ast (o-z)$ = 0)

for($w\leftarrow o$; $w\ge 0$; $w--$) $r\left[w\right]\leftarrow 0$; if($x=0$) $r\left[0\right]\leftarrow 1.0$ else $r\left[o\right]\leftarrow 1.0$; exit

EndIf

$r\left[z\right]\leftarrow z\ast \text{LOG}(z/o)+(o-z)\ast \text{LOG}(1-z/o)$

If($z<o/2$)

for($w\leftarrow z$; $w>0$; $w--$) $r\left[z\right]\leftarrow r\left[z\right]+\text{LOG}(o-z+w)-\text{LOG}\left(w\right)$

Else

for($w\leftarrow o-z$; $w>0$; $w--$) $r\left[z\right]\leftarrow r\left[z\right]+\text{LOG}(z+w)-\text{LOG}\left(w\right)$

EndIf

$r\left[z\right]\leftarrow \text{EXP}\left(r\right[z\left]\right)$

for($w\leftarrow z-1$; $w\ge 0$; $w--$) $r\left[w\right]\leftarrow r[w+1]\ast (o-z)\ast (w+1)/z/(o-w)$

for($w\leftarrow z+1$; $w\le m$; $w++$) $r\left[w\right]\leftarrow r[w-1]\ast x\ast (o-w+1)/w/(o-z)$

end procedure

Output:r    //$r\to f_{\mathrm{BS}}(0;x,m),\dots ,f_{\mathrm{BS}}(m;x,m)$ or $f_{\mathrm{BS}}(0;y,n),\dots ,f_{\mathrm{BS}}(n;y,n)$

Algorithm 2: Expression for ER, OR, RR (Equations (10)–(12)) as irreducible fraction.

Input:$x,m,y,n$    //x←x, m←m,y←y, n←n

procedure FR(a, b, &c, &d)

$c\leftarrow a$; $d\leftarrow b$;

If($d=0$) if($c=1$) $d=1$; $c=1$; exit EndIf

if($c=0$) $d\leftarrow 1$; if($d<2$) exit; if($a<0$) $a\leftarrow -a$

For(; $b>0$;)

$t\leftarrow b$; $b\leftarrow a\%b$; $a\leftarrow t$    //% is the modulo operator

EndFor    //a is the greatest common divisor

$c\leftarrow c/a$; $d\leftarrow d/a$    //$c/d$ is irreducible fraction

end procedure

procedure ER(x, y, m, n, &c, &d)    //excess risk, ER

FR(x*n - y*m, m*n, c, d)

end procedure   //$c/d\to$ ER (Equation (10))

procedure OR(x, y, m, n, &c, &d)    //odds ratio, OR

FR(x*(n-y), y*(m-x), c, d)

end procedure   //$c/d\to$ OR (Equation (11))

procedure RR(x, y, m, n, &c, &d)    //relative risk, RR

FR(x*n, y*m, c, d)

end procedure   //$c/d\to$ RR (Equation (12))

Output:$c,d$    //$\frac{c}{d}$ is a irreducible fraction

Algorithm 3: Bivariate binomial experiment, Equation (14).

Input:x, y, m, n, F    //x←x, y←y, m←m, n←n, $F\leftarrow$ ER, OR or RR

procedure B2F(x, y, m, n, F, &p, &g, &h)

PMF_B($x,m,r1$); PMF_B($y,n,r2$); $p\leftarrow \left[\right]$; $g\leftarrow \left[\right]$; $h\leftarrow \left[\right]$

For($u\leftarrow m;u\ge 0;u--$) For($v\leftarrow n;u\ge 0;v--$)

$p\left[\right]\leftarrow r1\left[u\right]\ast r2\left[v\right]$; $F(x,y,m,n,c,d)$; $g\left[\right]\leftarrow c$; $h\left[\right]\leftarrow d$

EndFor EndFor

end procedure

Output:$p,g,h$    //$p\left[\right]$ probabilities; $g\left[\right]$ numerators; $h\left[\right]$ denominators

Algorithm 4: PMF for bivariate binomial expressions (B2E), Equation (15).

Input:k, p, g, h    //k←k, p←p, g←g, h←h     procedure PMF_B2E(&k, &p, &g, &h)        procedure S0FR(l,r)         $i\leftarrow l$; $j\leftarrow r$; $t\leftarrow \text{DIV}(l+r,2)$         For(; $i>j$;)           for(; $g\left[i\right]/h\left[i\right]<g\left[t\right]/h\left[t\right]$; $i++$); for(; $g\left[t\right]/h\left[t\right]<g\left[j\right]/h\left[j\right]$; $j--$);           If($i\le j$)             $y\leftarrow g\left[i\right];g\left[i\right]\leftarrow g\left[j\right];g\left[j\right]\leftarrow y$; $y\leftarrow h\left[i\right];h\left[i\right]\leftarrow h\left[j\right];h\left[j\right]\leftarrow y$             $y\leftarrow p\left[i\right];p\left[i\right]\leftarrow p\left[j\right];p\left[j\right]\leftarrow y$; $i\leftarrow i+1;j\leftarrow j-1$           EndIf          EndFor       if($l<j$) $\text{S0FR}(l,j)$; if($i<r$) $\mathrm{S0FR}(i,r)$     end procedure    $\mathrm{S0FR}(0,k-1)$; $p1\leftarrow \left[\right]$; $g1\leftarrow \left[\right]$; $h1\leftarrow \left[\right]$; $k1\leftarrow 0$    For($i\leftarrow 0$; $i<k$;)       $s\leftarrow g\left[i\right]$       For($j\leftarrow i+1$; $j<k$;)          if($g\left[j\right]\ne g\left[i\right]$) $\text{BREAK}$; if($h\left[j\right]\ne h\left[i\right]$) $\text{BREAK}$; $s\leftarrow s+p\left[j\right]$       EndFor       $p1\left[\right]\leftarrow s$; $g1\left[\right]\leftarrow g\left[i\right]$; $h1\left[\right]\leftarrow h\left[i\right]$; $i\leftarrow j$; $k1\leftarrow k1+1$    EndFor    $k\leftarrow k1$; $p\leftarrow p1$; $g\leftarrow g1$; $h\leftarrow h1$ end procedure Output: $k,p,g,h$    //k→k, p→p, g→g, h→h

Algorithm 5: CIs in increasing coverage for B2E.

Input:$nu$, $de$, k, p, g, h    //$nu\leftarrow$nu, de←de, k←mn, p←p, g←g, h←h     function CC(&g, &h, k, i, j)     $a\leftarrow “[”$; $b\leftarrow “]”$     If($k>0$) if($i>0$) $i--$; $a\leftarrow “(”$ endif; if($j<k$) $j++$; $b\leftarrow “)”$ endif EndIf     $u\leftarrow \text{CONCATENATE}\left(g\right[i],“/”,h[i\left]\right)$; if($h\left[i\right]=1$) $u\leftarrow g\left[0\right]$; if($h\left[i\right]=0$) $u\leftarrow 0$     $v\leftarrow \text{CONCATENATE}\left(g\right[j],“/”,h[j\left]\right)$; if($h\left[j\right]=1$) $v\leftarrow g\left[0\right]$; if($h\left[j\right]=0$) $v\leftarrow 0$     $w\leftarrow \text{CONCATENATE}(a,u,“,”,v,b)$; $\text{RETURN}\left(w\right)$     end function     procedure CI_B2E(nu, de, &k, &p, &g, &h, &ci)     $i0\leftarrow -1$; for($i\leftarrow k-1$; $i\ge 0$; $i--$) if(($nu=g\left[i\right]$) AND ($de=h\left[i\right]$)) $i0\leftarrow i$     $it\leftarrow 0$; $i1\leftarrow i0$; $i2\leftarrow i0$; $iq\leftarrow 0$; $r\leftarrow p\left[i0\right]$; $q\leftarrow 1.0-p\left[i0\right]$     For(; ;)        if($r<0.5$) $q\leftarrow 1-r$ else $r\leftarrow 1-q$; $c0\leftarrow CC(g,h,0,i1,i2)$; $c1\leftarrow CC(g,h,k,i1,i2)$        $ci\left[it\right]\leftarrow \text{ARRAY}(c0,c1,r,q)$; $it\leftarrow it+1$; if(($i1=0$) AND ($i2=k$)) $\text{BREAK}$        if($i1=0$) $i2\leftarrow i2+1$; $r\leftarrow r+p\left[i2\right]$; $q\leftarrow q-p\left[i2\right]$; $\text{CONTINUE}$ endif        if($i2=k$) $i1\leftarrow i1-1$; $r\leftarrow r+p\left[i2\right]$; $q\leftarrow q-p\left[i2\right]$; $\text{CONTINUE}$ endif        If($p[i1-1]=p[i2+1]$)           $i1\leftarrow i1-1$; $i2\leftarrow i2+1$; $r\leftarrow r+p\left[i1\right]+p\left[i2\right]$; $q\leftarrow q-p\left[i1\right]-p\left[i2\right]$; $\text{CONTINUE}$        EndIf        if($p[i1-1]>p[i2+1]$) $i1\leftarrow i1-1$; $r\leftarrow r+p\left[i1\right]$; $q\leftarrow q-p\left[i1\right]$; $\text{CONTINUE}$ endif        if($p[i1-1]<p[i2+1]$) $i2\leftarrow i2+1$; $r\leftarrow r+p\left[i2\right]$; $q\leftarrow q-p\left[i2\right]$; $\text{CONTINUE}$ endif     EndFor     $k\leftarrow it$     end procedure Output: $ci$    //$ci\left[i\right]$ array (ci_min, ci_max, coverage, 1-coverage)

3. Background

Since ratio scale statistics are widely known, those shall be used for comparison. A long time ago [13], Carl Friedrich Gauss wrote the formula for the distribution of the population statistic known presently as “Normal”-Gauss distribution. Some readers may believe that Gauss was the first to address this issue, which is not quite true. With a reasonable confidence, one should say that Newton’s foundational work [14] has inspired many people. Drawing from this, Bernoulli [15] put forward a piece of proof for the expansion formula of the Newton binomial; thus, the Binomial distribution was invented ($f_{\mathrm{B}}$, Equation (1), where ${\Gamma }(o+1)=o!$).

$$\begin{array}{c}1={(p+1-p)}^n=\sum _{k=0}^n\left(\genfrac{}{}{0pt}{}{n}{k}\right)p^k{(1-p)}^{n-k}=\sum _{k=0}^n{f}_{\mathrm{B}}(k;n,p),\\ f_{\mathrm{B}}(k;n,p)=\frac{{\Gamma }(n+1)}{{\Gamma }(k+1){\Gamma }(n-k+1)}{p}^k{(1-p)}^{n-k}\end{array}$$

(1)

The first to look for asymptotic formulas to be used for the approximation of Binomial distribution was De Moivre [16]. Merely approximately 100 years later, Gauss [13] looked over the Gaussian function in an ampler manner ($f_{\mathrm{G}}$, Equation (2); where $a={\sigma }^{-1}{\left(2\pi \right)}^{-1/2},b=\mu ,c=\sigma$ for Normal distribution, $f_{\mathrm{N}}$ in Equation (3)).

$$f_{\mathrm{G}}(x;a,b,c)=a·\exp \left(-\frac{{(x-b)}^2}{2c^2}\right)$$

(2)

$$\begin{array}{c}{f}_{\mathrm{N}}(x;\mu ,\sigma )=\frac{1}{\sigma {\left(2\pi \right)}^{1/2}}\exp \left(-\frac{{(x-\mu )}^2}{2{\sigma }^2}\right),\\ {\int }_{-\infty }^{\infty }{f}_{\mathrm{N}}(x;\mu ,\sigma )dx=1\end{array}$$

(3)

Occurring in most of the natural phenomena, the Gauss distribution is used to model observed data collected on a ratio scale. To elaborate, it is a continuous distribution, which means that it starts from the assumption that there is an infinite diversity in the population of the data, and in between any two distinct values, one more (distinct value) still exists. Anyone should recognize that it is a very strong requirement, since most of the instrumentation is normally trained to collect the data with no more than four significant digits (leading to a universe build up from populations of only thousands of distinct individuals). Nevertheless, sampling (in which the limited instrumentation filters the collected data with a limited number of digits) generally operates on small-sized drawings (even smaller than thousands) from the population. Consequently, the population distribution is actually never reconstituted from any number of repeated samples. In other words, by sampling, Normal distribution is never resampled. In fact, through sampling one sees the population distribution differently (than the one which actually exists). In summary, this is the sampling distribution. For the Gauss distribution, the first to propose its sampling distribution was Gosset [17], and that is Student’s t distribution (Equation (4)).

$$f_{\mathrm{t}}(x;\nu )=\frac{{\Gamma }\left(\frac{\nu +1}{2}\right)}{{\left(\nu \pi \right)}^{1/2}{\Gamma }\left(\frac{\nu }{2}\right)}{\left(1+\frac{x^2}{\nu }\right)}^{-\frac{\nu +1}{2}}$$

(4)

In a universe overwhelmed by its infinity, Gauss distribution is also what a proportion will follow (see Equation (5); see Figure 2—Student’s t and Binomial against Normal).

$$f_{\mathrm{N}}(x;0,1)-f_t(x;\nu )\stackrel{\nu \to \infty }{\to }0,f_{\mathrm{N}}(x;np,\sqrt{np(1-p)})-f_B(x;n,p)\stackrel{n\to \infty }{\to }0$$

(5)

The same problem—asymptotic formulas for very large populations—was addressed in Physics and Chemistry as well, when important progress was being made to the Maxwell–Boltzmann [18,19,20], Bose–Einstein [21,22,23] and Fermi–Dirac [23,24,25] distributions (all discrete, all with “continuity correction”—to borrow the binomial related language—approximations intended to work only for large numbers [26])—but those were very lucky cases, since the population of atoms and molecules is indeed enormous (${\mathrm{N}}_{\mathrm{A}},\mathrm{about}6·{10}^{23}$ molecules exist in normal conditions in a relatively small amount of one mol of substance; ${\mathrm{N}}_{\mathrm{A}}$—Avogadro’s [27] number). Still, it is very useful to learn how physicists arrived to the continuous case. Thus, extracting relevant equations from [20], the following is a short route:

• For a closed system with indistinguishable particles and K energy states, the event to observe a certain configuration ($n_1,\dots ,n_K$) is a multinomial distribution ([28]; Equation (6)):

  $$f_{\mathrm{M}}(n_1,\dots ,n_K;p_1,\dots ,p_K)=\frac{{\prod }_{i=1}^K{p}_i^{n_i}}{{\prod }_{i=1}^K{n}_i!}\left(\sum _{i=1}^K{n}_i\right)!$$

  (6)
• A closed system with no chemical changes [29] is under two constrains, $N={\sum }_{i=1}^K{n}_i$ and $E={\sum }_{i=1}^K{e}_i{n}_i$; by maximizing the likelihood [30], according to Lagrange multiplier’s method [31]: $\psi (n_i+1)=\ln \left(p_i\right)-\alpha -\beta e_i$ (where $\psi$ is the digamma function);
• If $n_i\gg 1$ (see [32]) then $\psi (n_i+1)-\ln \left(n_i\right)\stackrel{n_i\to \infty }{\to }\gamma$ (Euler’s gamma constant [33]) and the discrete Maxwell–Boltzmann distribution of particles by energies is obtained ($p_i$ in Equation (7); where ${\sum }_{j=1}^K{e}^{-\beta e_j}$ is the partition function (introducing the entropy [34]).
• For an ideal gas, the value of a speed component (for instance x) takes (hypothetically) any value from $(-\infty ,\infty )$, so (${\Sigma }\to \int$) continuous Maxwell–Boltzmann distribution of particles by speed components is obtained, which is a Gaussian ($f_{\mathrm{GMB}}$ in Equation (7));
• Ultimately, for a speed s constituted from J components, the formula for continuous Maxwell–Boltzmann distribution becomes ($f_{MBJ}$ in Equation (7); Equation (29) in [20]):

  $$p_i=\frac{\frac{n_i}{N}{e}^{-\beta e_i}}{{\sum }_{j=1}^K{e}^{-\beta e_j}},f_{\mathrm{MBG}}\left(x\right)=a{e}^{-b{x}^2},f_{\mathrm{MBJ}}(s;a,J)=\frac{{\pi }^{J/2}}{2^{J-1}{\Gamma }(J/2)}{s}^{J-1}{a}^J{e}^{-\pi a^2{s}^2}$$

  (7)

However, since at least the sample sizes are finite (if not the populations too), one needs to address the sample size here as well.

To cite Alan Agresti—the (U.S.) Statistician of year 2003—who recommends “using confidence methods for proportions, for the difference between proportions, and for odd ratios in three situations: when a lower bound on a coverage probability (CP) is desired, when an actual CP near the nominal level is desired, and when teaching in a classroom.” Following [35], one could say “this is considered true with ’a certain level of confidence”.

Calling for the use of the exact method for the calculation of the CI at the difference of two binomial proportions, the authors of [36] made progress in a series of studies [37,38,39,40,41], but again, in most of the cases, the inverse problem is needed to find the probability associated with a certain difference or ratio of two proportions [42,43,44,45,46,47].

4. Material and Method

Independently of the fact that the difference of the proportions was expressed in percentages or absolute values, its limits are to be derived from the model of the research question in case. Taking into consideration the case of one proportion, there are three distinct cases (where $\alpha$ is the risk of being in error):

• Continuous asymptotic (normal) case, when CI is given by the inverse of the cumulative distribution function also called the percent-point function or quantile function: $\pm \sigma ·z\left(\alpha \right)$ (where ${\sigma }^2$ is the variance); for the occurrence of the proportion in a sample of size n, CI is $\pm \sigma ·z\left(\alpha \right)/\sqrt{n}$;
• Continuous non-asymptotic (student) case, proportion coming from a measurement through sampling of size m, when CI is given with Student t-value instead of the z-value: $\pm \sigma ·t(\alpha ,m)$; for the incidence in a sample of size n, CI is $\pm \sigma ·t(\alpha ,m)/\sqrt{n}$;
• Discrete case, when a proportion comes from a sample in a (no matter how large) population, one needs to know the value of the binomial variable (x) and the sample size (m) in order to express its CI; it always becomes suitable to use the binomial distributed CI in the calculation; one of the first proposed CI is in [48] ($InvCD{F}_{\beta }(1-\alpha /2,r+1,m-r)$, where $f_{\beta }(x;{\alpha }_1,{\alpha }_2)={\Gamma }({\alpha }_1+{\alpha }_2)·{\Gamma }{\left({\alpha }_1\right)}^{-1}·{\Gamma }{\left({\alpha }_2\right)}^{-1}·x^{{\alpha }_1-1}·{(1-x)}^{{\alpha }_2-1}$), and $InvCD{F}_{\beta }$ is its inverse—providing a CI at least as large as it is supposed to be and being most of the time too large; since the discreteness transforms the matter of integration (back) into a matter of summation and the calculation of each CI is a problem of combinatorics; there is no unique solution. One alternative is given in [9] and another one in [49], both of those alternatives have no analytical formula, only recipes (algorithms) leading to the solution.

Thus, when a difference of two proportions is done, it is actually a full factorial design of the cases given above, each of the two proportions possibly falling in one of the three categories given above.

• At least one continuous asymptotic (normal) case. When comparing two true values in a (very) large population (of an unknown size), then both proportions fall into the asymptotic continuous case; it may be suitable to use the z-value in the calculation (for studies involving z-value, z-statistics, z-test, see [50,51,52,53]);
• At least one continuous non-asymptotic (student) case. When comparing two true values obtained from a large number of repeated samplings (of known size) in a (very) large population, then both proportions fall into the continuous but not asymptotic case; it is much more suitable (than the call for the previous case) to use the Student t-value (for some typical cases, see [54,55,56,57]);
• 2 × 2 discrete asymptotic case. When comparing with a true value in a (very) large population (of an known or unknown size), it is a typical case of involving an asymptotic method (see the methods given in [58,59,60]);
• 2 × 2 discrete case. In most of the cases, the results compared are coming from two samples (of necessary known sizes), when actually a bidimensional binomial distribution is involved and the probability for the difference is to be calculated from it [36]. This case is detailed in the next paragraph and in the remainder of the paper.

Under different arrangements (see a different arrangement in [61]), a contingency table ([62,63]; a confusion matrix in [64]; an error matrix in [65]) collects observed frequencies at the contingency of two multinomial distributed variables. The simplest case, 2 × 2 contingency, is for two dichotomies (fourfold table in [63]; 2 × 2 table in [66]).

The 2 × 2 contingencies are used in at least two different contexts (see Figure 3).

Any repeated draw u from the population in respect to the sample S1 (see Equation (1); see Figure 3) has a probability of occurrence $f_{\mathrm{S}1}$, and $f_{\mathrm{S}2}$ for v in respect to S2 (Equation (8)).

$$f_{\mathrm{BS}}(w;z,o)=\frac{o!}{z!(o-z)!}{\left(\frac{z}{o}\right)}^w{\left(1-\frac{z}{o}\right)}^{o-w},f_{\mathrm{S}1}\leftarrow f_{\mathrm{BS}}(u;x,m),f_{\mathrm{S}2}=f_{\mathrm{BS}}(v;y,n)$$

(8)

A replica $(u/m,v/n)$ of the binomial experiment $(x/m,y/n)$ will be produced with a probability (which is the bivariate binomial, with a probability given by the product of the two binomial probabilities $f_{\mathrm{S}1}$ and $f_{\mathrm{S}2}$), Equation (9):

$$f_{\mathrm{B}2\mathrm{E}}(u,v;x,y,m,n)=\frac{m!}{u!(m-u)!}{\left(\frac{x}{m}\right)}^u{\left(1-\frac{x}{m}\right)}^{m-u}\frac{n!}{v!(n-v)!}{\left(\frac{y}{n}\right)}^v{\left(1-\frac{y}{n}\right)}^{n-v}$$

(9)

Equation (9) can be perceived by some as slightly intimidating, but it actually expresses a state of facts: any event involving u and v as possible outcomes has a conditional probability of the already observed (x and y) outcomes.

Excess risk (ER, Equation (10)), Odds ratio (OR, Equation (11)) and Relative risk (RR, Equation (12)) are algebraic expressions of the two proportions ($x/m$ and $y/n$) followed up in the study (Figure 3):

$$\begin{array}{c}ER=ER(x,y;m,n)=\frac{x}{m}-\frac{y}{n},\end{array}$$

(10)

$$\begin{array}{c}OR=OR(x,y;m,n)=\frac{x(n-y)}{y(m-x)},\end{array}$$

(11)

$$\begin{array}{c}RR=RR(x,y;m,n)=\frac{xn}{ym}\end{array}$$

(12)

Since the samples are independent draws from the population, with no further means to provide estimates of the characteristics of interest in the population (the $x/m$ and $y/n$ proportions, are of positive outcomes, and are the characteristics of interest), each of the two proportions ($x/m$ and $y/n$) comes from typically binomial distributed variables (outcome x having as domain $\{0,1,\dots ,m\}$ and outcome y having as domain $\{0,1,\dots ,n\}$) and, thus, each follows a binomial distribution while together they generate a bidimensional (or bivariate) binomial distribution.

There are essential differences between the concepts. The absolute risk reduction reflects the probability of getting a disease in the first place, for instance, while reporting only relative risk reduction, one would run into the risk that readers exaggerate the effectiveness of a treatment [67].

5. Results and Discussion

As given in the Research design schema (Figure 1), the solution to the confidence interval and/or p-value problem should always follow a series of steps.

Calculation of the CI for the discrete case is a problem of combinatorics. Therefore, a paramount issue is to have an efficient (in terms of the execution time) and accurate (in terms of the calculation errors) algorithm for the calculation of the probability from a bivariate binomial distribution (Equation (10)). Several variants were tried and several strategies were tested. Below are given some acquired tips:

• It is not efficient to directly evaluate Equation (10); the complexity of the evaluation is $\mathcal{O}\left(m^2{n}^2\right)$ and, for large values of the samples sizes (m and n), one will never see the result of the direct evaluation of Equation (10) (if for $m=n=100$ the evaluation takes 4 s, then for $m=n=$ 10,000, it takes 40,000 s, which is over 11 h);
• It is not advisable to blindly use a recursion formula either (see, for example, the recursion formulas below, Equation (13)); for instance, for m = n = 300, x = 15 and y = 30, by using a recursion formula starting from $f_{\mathrm{B}2\mathrm{E}}(m,n;x,y,m,n)$ and going to $f_{\mathrm{B}2\mathrm{E}}(0,0;x,y,m,n)$ by using a series of recursion formulas (such as the ones in Equation (13)), one will evaluate everything to 0 in a double (8 bytes, 64 bits) floating point precision for the simple reason that $f_{\mathrm{B}2\mathrm{E}}(m,n;x,y,m,n)$ is already lower ($4.9\dots ·{10}^{-691}$) than the lowest absolute value on the scale ($4.9\dots ·{10}^{-324}$) and is evaluated to 0 on any platform operating on IEEE 754 double-storage format;
• It is inefficient to use the general formula Equation (10) for the evaluation since for each u (or v) the same half-part $f_{\mathrm{BS}}(u;x,m)$ (or $f_{\mathrm{BS}}(v;y,n)$, respectively) of it is used m (or n, respectively) times, thus inducing an inflation of the calculations by m (or n, respectively) times; it is much more convenient to evaluate (and store) the $f_{\mathrm{BS}}(u;x,m)$ and $f_{\mathrm{BS}}(v;y,n)$ probabilities first, and later to simply multiply them ($f_{\mathrm{B}2\mathrm{E}}(u,v;x,y,m,n)\leftarrow f_{\mathrm{BS}}(u;x,m)·f_{\mathrm{BS}}(v;y,n)$).

A major observation is that $f_{\mathrm{B}2\mathrm{E}}(u,v;x,y,m,n)$ is at maximum when $u=x$ and $v=y$ and this maximum is still a subunitary number. Thus, independently of the magnitude of the sample sizes, the probability will be concentrated around the point of the coordinates $(x,y)$ located in the $\{0,1,\dots ,m\}\times \{0,1,\dots ,n\}$ grid. Therefore, an iterative evaluation (of $f_{\mathrm{BS}}(u;x,m)$ and $f_{\mathrm{BS}}(v;y,n)$ probabilities, see above) beginning from $(x,y)$ and iterating (descending) to $(0,0)$ and (ascending to) (m, n) will propagate through repeated multiplications and divisions the smallest amount of error in the most significant part of the probability (see Figure 4). Thus, both of the following recurrences are used (Equation (13)):

$$\begin{array}{c}\hfill f_{\mathrm{BS}}(u;x,m)=f_{\mathrm{BS}}(u-1;x,m)\frac{x(m-u+1)}{u(m-x)},f_{\mathrm{BS}}(u;x,m)=f_{\mathrm{BS}}(u+1;x,m)\frac{(m-x)(u+1)}{x(m-u)}\\ \hfill f_{\mathrm{BS}}(v;y,n)=f_{\mathrm{BS}}(v-1;y,n)\frac{y(n-v+1)}{v(n-y)},f_{\mathrm{BS}}(v;y,n)=f_{\mathrm{BS}}(v+1;y,n)\frac{(n-y)(v+1)}{y(n-v)}\end{array}$$

(13)

5.1. The Expressions of Two Binomials and Their CIs

Following up the notation from Figure 3, x defines the number of successes occurring in a number of m successive draws from a population, while y defines the number of successes occurring in a number of n successive draws. The bridge between the two is the population, supposedly being the same, each of the two experiments following a different characteristic of interest. Therefore, each sample comes from a binomial experiment, while together they form a bidimensional or bivariate binomial experiment.

A draw in the sample S1 (or in any of its replica; following up one characteristic in the population) is independent from a draw in the sample S2 (or in any of its replica; following up (the) other characteristic in the population). Thus, if $(u,v,m,n)$ is a replica of the $(x,y,m,n)$ experiment, then it is safely decomposed in two binomial ones, $(u,m)$ as replica of $(x,m)$ (having $m+1$ possible outcomes, from $u=0$ to $u=m$), and $(v,n)$ as replica of $(y,n)$ (having $n+1$ possible outcomes, from $u=0$ to $u=m$), each event having associated a probability (see

Table 1. Probability mass functions for the replica ($u,v,m,n$) of the 2 × 2 contingency ($x,y,m,n$).

Events	Probability Mass Function
($u,m$) replica of ($x,m$)	$\{f_{\mathrm{BS}}(0;x,m)$, $f_{\mathrm{BS}}(1;x,m)$, …, $f_{\mathrm{BS}}(u;x,m)$, …, $f_{\mathrm{BS}}(m;x,m)\}$
($v,n$) replica of ($y,n$)	$\{f_{\mathrm{BS}}(0;y,n)$, $f_{\mathrm{BS}}(1;y,n)$, …, $f_{\mathrm{BS}}(v;y,n)$, …, $f_{\mathrm{BS}}(n;y,n)\}$

For the expressions of $f_{\mathrm{BS}}(u;x,m)$ and $f_{\mathrm{BS}}(v;y,n)$, see Equation (8). The probability for $(u,v,m,n)$ as replica of $(x,y,m,n)$ is, again, $f_{\mathrm{B}2\mathrm{E}}(u,v;x,y,m,n)$ from Equation (9).

).

One should know that following the design of Table 1, the CI will be provided for the sampled (the one collected in the 2 × 2 contingency, see Figure 3) with no assumption of the expected value in the population, or, more importantly, if the population is finite or infinite (usually an assumption leading to ’asymptotic CIs’). Turning back to the Gauss vs. Student challenge from the beginning of this paper (Equation (3) vs. Equation (4)), one could say that what is proposed in this case is a Student-type correction (for a finite repeated resampling) against an infinite (size) Gauss sampling.

An example will be provided and continued below. For $x=2$, $y=1$, $m=6$ and $n=4$ as a particular case of a bivariate binomial experiment, the probability mass functions are given in

Table 2. Example: probability mass functions for x = 2, $y=1$, m = 6 and n = 4 bivariate binomial experiment case.

${\mathit{f}}_{\mathbf{BS}}$	0	1	2	3	4	5	6
$f_{\mathrm{BS}}(u;2,6)$	$8.78·{10}^{-2}$	$2.63·{10}^{-1}$	$3.29·{10}^{-1}$	$2.19·{10}^{-1}$	$8.23·{10}^{-2}$	$1.65·{10}^{-2}$	$1.37·{10}^{-3}$
$f_{\mathrm{BS}}(v;1,4)$	$3.16·{10}^{-1}$	$4.22·{10}^{-1}$	$2.11·{10}^{-1}$	$4.69·{10}^{-2}$	$3.91·{10}^{-3}$

The numeric values are given with three significant figures. However, the calculations and the actual values must be done with machine-like precision. BS($u;2,6$) has seven possible outcomes ($u\in \{0,1,2,3,4,5,6\}$) while BS($v;1,4$) has five ($v\in \{0,1,2,3,4\}$).

.

Probability mass functions for expressions of two proportions. If $\mathrm{B}2\mathrm{E}(x,y,m,n)$ is the bivariate binomial experiment, then it has $(m+1)\times (n+1)$ possible outcomes (not necessarily seen distinctly through an evaluation function as the ones given in Equations (10)–(12) but necessary), each controlled by a probability given by $f_{\mathrm{B}2\mathrm{E}}(u,v;x,y,m,n)$. A list of $(m+1)\times (n+1)$ paired values should be collected at this stage. By using a simplified notation: $f_{1,2}(u,v)=f_{\mathrm{B}2\mathrm{E}}(u,v;x,y,m,n)$ and $h_{1,2}(x,y)=F2(x,y,m,n)$, where $F2(x,y,m,n)$ is a generic notation for any expression involving two binomial proportions (including the ones given in Equations (10)–(12)), the result of the pairing is as given below (Equation (14)):

$$\begin{array}{c}\left\{\right(h_{1,2}(0,0),f_{1,2}(0,0),\dots ,(h_{1,2}(0,v),f_{1,2}(0,v)),\dots ,(h_{1,2}(0,n),f_{1,2}(0,n)),\\ \dots \\ (h_{1,2}(u,0),f_{1,2}(u,0)),\dots .,(h_{1,2}(u,v),f_{1,2}(u,v)),\dots ,(h_{1,2}(u,n),f_{1,2}(u,n)),\\ \dots \\ (h_{1,2}(m,0),f_{1,2}(m,0)),\dots ,(h_{1,2}(m,v),f_{1,2}(m,v)),\dots ,(h_{1,2}(m,n),f_{1,2}(m,n))\}\end{array}$$

(14)

Following up the example in Table 2, the Appendix

Table A1. Pairs of possible ERs and their associated probabilities ($h_{1,2},f_{1,2}$) for the data in Table 2.

u vs. v	0	1	2	3	4
0	($\frac{0}{6}-\frac{0}{4}$, $2.78·{10}^{-2}$)	($\frac{0}{6}-\frac{1}{4}$, $3.70·{10}^{-2}$)	($\frac{0}{6}-\frac{2}{4}$, $1.85·{10}^{-2}$)	($\frac{0}{6}-\frac{3}{4}$, $4.12·{10}^{-3}$)	($\frac{0}{6}-\frac{4}{4}$, $3.43·{10}^{-4}$)
1	($\frac{1}{6}-\frac{0}{4}$, $8.33·{10}^{-2}$)	($\frac{1}{6}-\frac{1}{4}$, $1.11·{10}^{-1}$)	($\frac{1}{6}-\frac{2}{4}$, $5.56·{10}^{-2}$)	($\frac{1}{6}-\frac{3}{4}$, $1.23·{10}^{-2}$)	($\frac{1}{6}-\frac{4}{4}$, $1.03·{10}^{-3}$)
2	($\frac{2}{6}-\frac{0}{4}$, $1.04·{10}^{-1}$)	($\frac{2}{6}-\frac{1}{4}$, $1.39·{10}^{-1}$)	($\frac{2}{6}-\frac{2}{4}$, $6.94·{10}^{-2}$)	($\frac{2}{6}-\frac{3}{4}$, $1.54·{10}^{-2}$)	($\frac{2}{6}-\frac{4}{4}$, $1.29·{10}^{-3}$)
3	($\frac{3}{6}-\frac{0}{4}$, $6.94·{10}^{-2}$)	($\frac{3}{6}-\frac{1}{4}$, $9.26·{10}^{-2}$)	($\frac{3}{6}-\frac{2}{4}$, $4.63·{10}^{-2}$)	($\frac{3}{6}-\frac{3}{4}$, $1.03·{10}^{-2}$)	($\frac{3}{6}-\frac{4}{4}$, $8.57·{10}^{-4}$)
4	($\frac{4}{6}-\frac{0}{4}$, $2.60·{10}^{-2}$)	($\frac{4}{6}-\frac{1}{4}$, $3.47·{10}^{-2}$)	($\frac{4}{6}-\frac{2}{4}$, $1.74·{10}^{-2}$)	($\frac{4}{6}-\frac{3}{4}$, $3.86·{10}^{-3}$)	($\frac{4}{6}-\frac{4}{4}$, $3.22·{10}^{-4}$)
5	($\frac{5}{6}-\frac{0}{4}$, $5.21·{10}^{-3}$)	($\frac{5}{6}-\frac{1}{4}$, $6.94·{10}^{-3}$)	($\frac{5}{6}-\frac{2}{4}$, $3.47·{10}^{-3}$)	($\frac{5}{6}-\frac{3}{4}$, $7.72·{10}^{-4}$)	($\frac{5}{6}-\frac{4}{4}$, $6.43·{10}^{-5}$)
6	($\frac{6}{6}-\frac{0}{4}$, $4.34·{10}^{-4}$)	($\frac{6}{6}-\frac{1}{4}$, $5.79·{10}^{-4}$)	($\frac{6}{6}-\frac{2}{4}$, $2.89·{10}^{-4}$)	($\frac{6}{6}-\frac{3}{4}$, $6.43·{10}^{-5}$)	($\frac{6}{6}-\frac{4}{4}$, $5.36·{10}^{-6}$)

gives the ER and their associated probabilities while

Table A2. Pairs of possible ORs and their associated probabilities ($h_{1,2},f_{1,2}$) for the data in Table 2.

u vs. v	0	1	2	3	4
0	($\frac{0\ast 4}{0\ast 6}$, $2.78·{10}^{-2}$)	($\frac{0\ast 3}{1\ast 6}$, $3.70·{10}^{-2}$)	($\frac{0\ast 2}{2\ast 6}$, $1.85·{10}^{-2}$)	($\frac{0\ast 1}{3\ast 6}$, $4.12·{10}^{-3}$)	($\frac{0\ast 0}{4\ast 6}$, $3.43·{10}^{-4}$)
1	($\frac{1\ast 4}{0\ast 5}$, $8.33·{10}^{-2}$)	($\frac{1\ast 3}{1\ast 5}$, $1.11·{10}^{-1}$)	($\frac{1\ast 2}{2\ast 5}$, $5.56·{10}^{-2}$)	($\frac{1\ast 1}{3\ast 5}$, $1.23·{10}^{-2}$)	($\frac{1\ast 0}{4\ast 5}$, $1.03·{10}^{-3}$)
2	($\frac{2\ast 4}{0\ast 4}$, $1.04·{10}^{-1}$)	($\frac{2\ast 3}{1\ast 4}$, $1.39·{10}^{-1}$)	($\frac{2\ast 2}{2\ast 4}$, $6.94·{10}^{-2}$)	($\frac{2\ast 1}{3\ast 4}$, $1.54·{10}^{-2}$)	($\frac{2\ast 0}{4\ast 4}$, $1.29·{10}^{-3}$)
3	($\frac{3\ast 4}{0\ast 3}$, $6.94·{10}^{-2}$)	($\frac{3\ast 3}{1\ast 3}$, $9.26·{10}^{-2}$)	($\frac{3\ast 2}{2\ast 3}$, $4.63·{10}^{-2}$)	($\frac{3\ast 1}{3\ast 3}$, $1.03·{10}^{-2}$)	($\frac{3\ast 0}{4\ast 3}$, $8.57·{10}^{-4}$)
4	($\frac{4\ast 4}{0\ast 2}$, $2.60·{10}^{-2}$)	($\frac{4\ast 3}{1\ast 2}$, $3.47·{10}^{-2}$)	($\frac{4\ast 2}{2\ast 2}$, $1.74·{10}^{-2}$)	($\frac{4\ast 1}{3\ast 2}$, $3.86·{10}^{-3}$)	($\frac{4\ast 0}{4\ast 2}$, $3.22·{10}^{-4}$)
5	($\frac{5\ast 4}{0\ast 1}$, $5.21·{10}^{-3}$)	($\frac{5\ast 3}{1\ast 1}$, $6.94·{10}^{-3}$)	($\frac{5\ast 2}{2\ast 1}$, $3.47·{10}^{-3}$)	($\frac{5\ast 1}{3\ast 1}$, $7.72·{10}^{-4}$)	($\frac{5\ast 0}{4\ast 1}$, $6.43·{10}^{-5}$)
6	($\frac{6\ast 4}{0\ast 0}$, $4.34·{10}^{-4}$)	($\frac{6\ast 3}{1\ast 0}$, $5.79·{10}^{-4}$)	($\frac{6\ast 2}{2\ast 0}$, $2.89·{10}^{-4}$)	($\frac{6\ast 1}{3\ast 0}$, $6.43·{10}^{-5}$)	($\frac{6\ast 0}{4\ast 0}$, $5.36·{10}^{-6}$)

and

Table A3. Pairs of possible RRs and their associated probabilities ($h_{1,2},f_{1,2}$) for the data in Table 2.

u vs. v	0	1	2	3	4
0	($\frac{0\ast 4}{0\ast 6}$, $2.78·{10}^{-2}$)	($\frac{0\ast 4}{1\ast 6}$, $3.70·{10}^{-2}$)	($\frac{0\ast 4}{2\ast 6}$, $1.85·{10}^{-2}$)	($\frac{0\ast 4}{3\ast 6}$, $4.12·{10}^{-3}$)	($\frac{0\ast 4}{4\ast 6}$, $3.43·{10}^{-4}$)
1	($\frac{1\ast 4}{0\ast 6}$, $8.33·{10}^{-2}$)	($\frac{1\ast 4}{1\ast 6}$, $1.11·{10}^{-1}$)	($\frac{1\ast 4}{2\ast 6}$, $5.56·{10}^{-2}$)	($\frac{1\ast 4}{3\ast 6}$, $1.23·{10}^{-2}$)	($\frac{1\ast 4}{4\ast 6}$, $1.03·{10}^{-3}$)
2	($\frac{2\ast 4}{0\ast 6}$, $1.04·{10}^{-1}$)	($\frac{2\ast 4}{1\ast 6}$, $1.39·{10}^{-1}$)	($\frac{2\ast 4}{2\ast 6}$, $6.94·{10}^{-2}$)	($\frac{2\ast 4}{3\ast 6}$, $1.54·{10}^{-2}$)	($\frac{2\ast 4}{4\ast 6}$, $1.29·{10}^{-3}$)
3	($\frac{3\ast 4}{0\ast 6}$, $6.94·{10}^{-2}$)	($\frac{3\ast 4}{1\ast 6}$, $9.26·{10}^{-2}$)	($\frac{3\ast 4}{2\ast 6}$, $4.63·{10}^{-2}$)	($\frac{3\ast 4}{3\ast 6}$, $1.03·{10}^{-2}$)	($\frac{3\ast 4}{4\ast 6}$, $8.57·{10}^{-4}$)
4	($\frac{4\ast 4}{0\ast 6}$, $2.60·{10}^{-2}$)	($\frac{4\ast 4}{1\ast 6}$, $3.47·{10}^{-2}$)	($\frac{4\ast 4}{2\ast 6}$, $1.74·{10}^{-2}$)	($\frac{4\ast 4}{3\ast 6}$, $3.86·{10}^{-3}$)	($\frac{4\ast 4}{4\ast 6}$, $3.22·{10}^{-4}$)
5	($\frac{5\ast 4}{0\ast 6}$, $5.21·{10}^{-3}$)	($\frac{5\ast 4}{1\ast 6}$, $6.94·{10}^{-3}$)	($\frac{5\ast 4}{2\ast 6}$, $3.47·{10}^{-3}$)	($\frac{5\ast 4}{3\ast 6}$, $7.72·{10}^{-4}$)	($\frac{5\ast 4}{4\ast 6}$, $6.43·{10}^{-5}$)
6	($\frac{6\ast 4}{0\ast 6}$, $4.34·{10}^{-4}$)	($\frac{6\ast 4}{1\ast 6}$, $5.79·{10}^{-4}$)	($\frac{6\ast 4}{2\ast 6}$, $2.89·{10}^{-4}$)	($\frac{6\ast 4}{3\ast 6}$, $6.43·{10}^{-5}$)	($\frac{6\ast 4}{4\ast 6}$, $5.36·{10}^{-6}$)

do the same for OR and RR.

It may be stated (as a Theorem) that, in the previously given list, for any equation from (10) to (12), there is always a duplicate but the proof is immediate if one checks that $ER(0,0,m,n)$ = 0 = $ER(m,n,m,n)$ (for ER), $OR(0,v,m,n)$ = 0 = $OR(u,n,m,n)$ (for OR, for any $u,v\ne 0$), and $RR(0,v,m,n)$ = 0 (for RR, for any $v\ne 0$). Undoubtedly, the number of the duplicates is favored by the common factors in the values of the sample sizes. It can be effortlessly proved that in the same list above there are at least $2·\max (m,n)$ values (and the minimum is reached when $m=n$). This fact supports the necessity of grouping (of the above given pairs) by the value of calculated expression ($h_{1,2}$ function above) when the probabilities ($f_{1,2}$ function above) are summed. To ease up the grouping process, but also as is needed later, the resulted list of pairs is sorted (ascending) by the values of the calculated expression. At this point, it is no longer relevant which values of the freshly sampled binomial variates (u and v) produced the outcomes; thus, the values in this new list are labeled by their position (from 1 to k, k being the number of $h_{1,2}(u,v)$ distinct values), Equation (15):

$$\{(h_1,f_1),\dots ,(h_k,f_k)\}$$

(15)

The grouping and sorting for the example case is illustrated in

Table A4. Sorted and grouped ERs ($h_j$) and their associated probabilities ($f_j$) for the data in Table A1.

j	${\mathit{h}}_{\mathit{j}}$	${\mathit{f}}_{\mathit{j}}$	Group j	j	${\mathit{h}}_{\mathit{j}}$	${\mathit{f}}_{\mathit{j}}$	Group j	j	${\mathit{h}}_{\mathit{j}}$	${\mathit{f}}_{\mathit{j}}$	Group j
1	$\frac{-1}{1}$	$3.43·{10}^{-4}$	{$\frac{0}{6}-\frac{4}{4}$}	9	$\frac{-1}{4}$	$4.73·{10}^{-2}$	{$\frac{3}{6}-\frac{3}{4}$, $\frac{0}{6}-\frac{1}{4}$}	17	$\frac{5}{12}$	$3.47·{10}^{-2}$	{$\frac{4}{6}-\frac{1}{4}$}
2	$\frac{-5}{6}$	$1.03·{10}^{-3}$	{$\frac{1}{6}-\frac{4}{4}$}	10	$\frac{-1}{6}$	$6.95·{10}^{-2}$	{$\frac{5}{6}-\frac{4}{4}$, $\frac{2}{6}-\frac{2}{4}$}	18	$\frac{1}{2}$	$6.97·{10}^{-2}$	{$\frac{6}{6}-\frac{2}{4}$, $\frac{3}{6}-\frac{0}{4}$}
3	$\frac{-3}{4}$	$4.12·{10}^{-3}$	{$\frac{0}{6}-\frac{3}{4}$}	11	$\frac{-1}{12}$	$1.15·{10}^{-1}$	{$\frac{4}{6}-\frac{3}{4}$, $\frac{1}{6}-\frac{1}{4}$}	19	$\frac{7}{12}$	$6.94·{10}^{-3}$	{$\frac{5}{6}-\frac{1}{4}$}
4	$\frac{-2}{3}$	$1.29·{10}^{-3}$	{$\frac{2}{6}-\frac{4}{4}$}	12	$\frac{0}{1}$	$7.41·{10}^{-2}$	{$\frac{6}{6}-\frac{4}{4}$, $\frac{0}{6}-\frac{0}{4}$, $\frac{3}{6}-\frac{2}{4}$}	20	$\frac{2}{3}$	$2.60·{10}^{-2}$	{$\frac{4}{6}-\frac{0}{4}$}
5	$\frac{-7}{12}$	$1.23·{10}^{-2}$	{$\frac{1}{6}-\frac{3}{4}$}	13	$\frac{1}{12}$	$1.40·{10}^{-1}$	{$\frac{5}{6}-\frac{3}{4}$, $\frac{2}{6}-\frac{1}{4}$}	21	$\frac{3}{4}$	$5.79·{10}^{-4}$	{$\frac{6}{6}-\frac{1}{4}$}
6	$\frac{-1}{2}$	$1.94·{10}^{-2}$	{$\frac{3}{6}-\frac{4}{4}$, $\frac{0}{6}-\frac{2}{4}$}	14	$\frac{1}{6}$	$1.01·{10}^{-1}$	{$\frac{4}{6}-\frac{2}{4}$, $\frac{1}{6}-\frac{0}{4}$}	22	$\frac{5}{6}$	$5.21·{10}^{-3}$	{$\frac{5}{6}-\frac{0}{4}$}
7	$\frac{-5}{12}$	$1.54·{10}^{-2}$	{$\frac{2}{6}-\frac{3}{4}$}	15	$\frac{1}{4}$	$9.27·{10}^{-2}$	{$\frac{6}{6}-\frac{3}{4}$, $\frac{3}{6}-\frac{1}{4}$}	23	$\frac{1}{1}$	$4.34·{10}^{-4}$	{$\frac{6}{6}-\frac{0}{4}$}
8	$\frac{-1}{3}$	$5.59·{10}^{-2}$	{$\frac{4}{6}-\frac{4}{4}$, $\frac{1}{6}-\frac{2}{4}$}	16	$\frac{1}{3}$	$1.08·{10}^{-1}$	{$\frac{5}{6}-\frac{2}{4}$, $\frac{2}{6}-\frac{0}{4}$}

for ER, in

Table A5. Sorted and grouped ORs ($h_j$) and their associated probabilities ($f_j$) for the data in Table A2.

j	${\mathit{h}}_{\mathit{j}}$	${\mathit{f}}_{\mathit{j}}$	Group j	j	${\mathit{h}}_{\mathit{j}}$	${\mathit{f}}_{\mathit{j}}$	Group j	j	${\mathit{h}}_{\mathit{j}}$	${\mathit{f}}_{\mathit{j}}$	Group j
1	$\frac{0}{1}$	$6.36·{10}^{-2}$	Group 0	7	$\frac{3}{5}$	$1.11·{10}^{-1}$	{$\frac{1\ast 3}{1\ast 5}$}	13	$\frac{3}{1}$	$9.26·{10}^{-2}$	{$\frac{3\ast 3}{1\ast 3}$}
2	$\frac{1}{15}$	$1.23·{10}^{-2}$	{$\frac{1\ast 1}{3\ast 5}$}	8	$\frac{2}{3}$	$3.86·{10}^{-3}$	{$\frac{4\ast 1}{3\ast 2}$}	14	$\frac{5}{1}$	$3.47·{10}^{-3}$	{$\frac{5\ast 2}{2\ast 1}$}
3	$\frac{1}{6}$	$1.54·{10}^{-2}$	{$\frac{2\ast 1}{3\ast 4}$}	9	$\frac{1}{1}$	$7.41·{10}^{-2}$	{$\frac{6\ast 0}{4\ast 0}$, $\frac{0\ast 4}{0\ast 6}$, $\frac{3\ast 2}{2\ast 3}$}	15	$\frac{6}{1}$	$3.47·{10}^{-2}$	{$\frac{4\ast 3}{1\ast 2}$}
4	$\frac{1}{5}$	$5.56·{10}^{-2}$	{$\frac{1\ast 2}{2\ast 5}$}	10	$\frac{3}{2}$	$1.39·{10}^{-1}$	{$\frac{2\ast 3}{1\ast 4}$}	16	$\frac{15}{1}$	$6.94·{10}^{-3}$	{$\frac{5\ast 3}{1\ast 1}$}
5	$\frac{1}{3}$	$1.03·{10}^{-2}$	{$\frac{3\ast 1}{3\ast 3}$}	11	$\frac{5}{3}$	$7.72·{10}^{-4}$	{$\frac{5\ast 1}{3\ast 1}$}	17	$\frac{1}{0}$	$2.90·{10}^{-1}$	Group 17
6	$\frac{1}{2}$	$6.94·{10}^{-2}$	{$\frac{2\ast 2}{2\ast 4}$}	12	$\frac{2}{1}$	$1.74·{10}^{-2}$	{$\frac{4\ast 2}{2\ast 2}$}

Group 0 is {$\frac{5\ast 0}{4\ast 1}$, $\frac{4\ast 0}{4\ast 2}$, $\frac{0\ast 0}{4\ast 6}$, $\frac{3\ast 0}{4\ast 3}$, $\frac{1\ast 0}{4\ast 5}$, $\frac{2\ast 0}{4\ast 4}$, $\frac{0\ast 1}{3\ast 6}$, $\frac{0\ast 2}{2\ast 6}$, $\frac{0\ast 3}{1\ast 6}$}. Group 17 is {$\frac{6\ast 1}{3\ast 0}$, $\frac{6\ast 2}{2\ast 0}$, $\frac{6\ast 4}{0\ast 0}$, $\frac{6\ast 3}{1\ast 0}$, $\frac{5\ast 4}{0\ast 1}$, $\frac{4\ast 4}{0\ast 2}$, $\frac{3\ast 4}{0\ast 3}$, $\frac{1\ast 4}{0\ast 5}$, $\frac{2\ast 4}{0\ast 4}$}.

for OR and

Table A6. Sorted and grouped relative risks ($h_j$) and their associated probabilities ($f_j$) for the data in Table A3.

j	${\mathit{h}}_{\mathit{j}}$	${\mathit{f}}_{\mathit{j}}$	Group j	j	${\mathit{h}}_{\mathit{j}}$	${\mathit{f}}_{\mathit{j}}$	Group j	j	${\mathit{h}}_{\mathit{j}}$	${\mathit{f}}_{\mathit{j}}$	Group j
1	$\frac{0}{1}$	$6.00·{10}^{-2}$	{$\frac{0\ast 4}{4\ast 6}$, $\frac{0\ast 4}{3\ast 6}$, $\frac{0\ast 4}{2\ast 6}$, $\frac{0\ast 4}{1\ast 6}$}	7	$\frac{2}{3}$	$1.91·{10}^{-1}$	{$\frac{4\ast 4}{4\ast 6}$, $\frac{3\ast 4}{3\ast 6}$, $\frac{2\ast 4}{2\ast 6}$, $\frac{1\ast 4}{1\ast 6}$}	13	$\frac{5}{3}$	$3.47·{10}^{-3}$	{$\frac{5\ast 4}{2\ast 6}$}
2	$\frac{1}{6}$	$1.03·{10}^{-3}$	{$\frac{1\ast 4}{4\ast 6}$}	8	$\frac{5}{6}$	$6.43·{10}^{-5}$	{$\frac{5\ast 4}{4\ast 6}$}	14	$\frac{2}{1}$	$9.29·{10}^{-2}$	{$\frac{6\ast 4}{2\ast 6}$, $\frac{3\ast 4}{1\ast 6}$}
3	$\frac{2}{9}$	$1.23·{10}^{-2}$	{$\frac{1\ast 4}{3\ast 6}$}	9	$\frac{8}{9}$	$3.86·{10}^{-3}$	{$\frac{4\ast 4}{3\ast 6}$}	15	$\frac{8}{3}$	$3.47·{10}^{-2}$	{$\frac{4\ast 4}{1\ast 6}$}
4	$\frac{1}{3}$	$5.68·{10}^{-2}$	{$\frac{2\ast 4}{4\ast 6}$, $\frac{1\ast 4}{2\ast 6}$}	10	$\frac{1}{1}$	$7.41·{10}^{-2}$	{$\frac{6\ast 4}{4\ast 6}$, $\frac{0\ast 4}{0\ast 6}$, $\frac{3\ast 4}{2\ast 6}$}	16	$\frac{10}{3}$	$6.94·{10}^{-3}$	{$\frac{5\ast 4}{1\ast 6}$}
5	$\frac{4}{9}$	$1.54·{10}^{-2}$	{$\frac{2\ast 4}{3\ast 6}$}	11	$\frac{10}{9}$	$7.72·{10}^{-4}$	{$\frac{5\ast 4}{3\ast 6}$}	17	$\frac{4}{1}$	$5.79·{10}^{-4}$	{$\frac{6\ast 4}{1\ast 6}$}
6	$\frac{1}{2}$	$8.57·{10}^{-4}$	{$\frac{3\ast 4}{4\ast 6}$}	12	$\frac{4}{3}$	$1.56·{10}^{-1}$	{$\frac{6\ast 4}{3\ast 6}$, $\frac{4\ast 4}{2\ast 6}$, $\frac{2\ast 4}{1\ast 6}$}	18	$\frac{1}{0}$	$2.89·{10}^{-1}$	{$\frac{6\ast 4}{0\ast 6}$, $\frac{5\ast 4}{0\ast 6}$, $\frac{4\ast 4}{0\ast 6}$, $\frac{3\ast 4}{0\ast 6}$, $\frac{1\ast 4}{0\ast 6}$, $\frac{2\ast 4}{0\ast 6}$}

for RR (to be found in the Appendix A).

Cumulative distribution functions and CIs. Equation (15) (illustrated for the selected case in Table A4, Table A5 and Table A6) is the cornerstone for any CI evaluation associated with any expression of two binomial proportions (including Equations (10)–(12)). Thus,

• If one wants a CI of $\alpha$ (say, 5%) risk of being in error, then it may start from the pair cumulating the probability for $u=x$ and $v=y$ (it will always be a unique entry for it in $(h_i,f_i)$, $i=1,\dots ,k$ list, the one for which $h_i$ is the irreducible fraction of the expression subjected to analysis). Let’s label that index with $i_0$ (for instance $i_0$ = 13 in Table A4 since ER($2,6,1,4$) = $\frac{1}{12}$, $i_0$ = 10 in Table A5 since OR($2,1,6,4$) = $\frac{3}{2}$, and $i_0$ = 12 in Table A6 since RR($2,1,6,4$) = $\frac{4}{3}$). Let us use a variable to cumulate the probability, $p_{i_1,i_2}$ ($p_{i_1,i_2}\leftarrow 0$). The open interval ($h_{i_0-1}$, $h_{i_0+1}$) contains exactly one possibility (when $u=x$ and $v=y$) and has the probability to appear $f_{i_0}$; thus, $p_{i_1,i_2}\leftarrow p_{i_1,i_2}+f_{i_0}$. By using two pointers ($i_1\leftarrow i_0$ and $i_2\leftarrow i_0$) to move down ($i_1\leftarrow i_1-1$) and up ($i_2\leftarrow i_2+1$) the CI ($h_{i_1}$, $h_{i_2}$) is expanded (when $i_1\leftarrow i_1-1$ also $p_{i_1,i_2}\leftarrow p_{i_1,i_2}+f_{i_1}$, and when $i_2\leftarrow i_2+1$ also $p_{i_1,i_2}\leftarrow p_{i_1,i_2}+f_{i_2}$) until (and no further than) $1-p_{i_1,i_2}<\alpha$. Regarding the data given in Table A4, Table A5 and Table A6 as example, the $1-\alpha$ CI is to be obtained (in each case) by summing the $f_i$ values for a block delimited by two boundaries; all the work is exemplified in the (next) series of three tables (

  Table 3. CIs in increasing coverage for ER (data from Table A4).

  i	MinCI	MaxCI	pCI	1 − pCI	i	MinCI	MaxCI	pCI	1 − pCI
  0	$[\frac{1}{12},\frac{1}{12}]$	$(0,\frac{1}{6})$	0.1397	0.8603	12	$[-\frac{1}{2},\frac{1}{2}]$	$(-\frac{7}{12},\frac{7}{12})$	0.9417	0.0583
  1	$[\frac{1}{12},\frac{1}{6}]$	$(0,\frac{1}{4})$	0.2404	0.7596	13	$[-\frac{7}{12},\frac{1}{2}]$	$(-\frac{2}{3},\frac{7}{12})$	0.9540	0.0460
  2	$[\frac{1}{12},\frac{1}{4}]$	$(0,\frac{1}{3})$	0.3330	0.6670	14	$[-\frac{7}{12},\frac{7}{12}]$	$(-\frac{2}{3},\frac{2}{3})$	0.9610	0.0390
  3	$[\frac{1}{12},\frac{1}{3}]$	$(0,\frac{5}{12})$	0.4407	0.5593	15	$[-\frac{7}{12},\frac{2}{3}]$	$(-\frac{2}{3},\frac{3}{4})$	0.9870	0.0130
  4	$[0,\frac{1}{3}]$	$(-\frac{1}{12},\frac{5}{12})$	0.5147	0.4853	16	$[-\frac{2}{3},\frac{2}{3}]$	$(-\frac{3}{4},\frac{3}{4})$	0.9883	0.0117
  5	$[-\frac{1}{12},\frac{1}{3}]$	$(-\frac{1}{6},\frac{5}{12})$	0.6297	0.3703	17	$[-\frac{3}{4},\frac{2}{3}]$	$(-\frac{5}{6},\frac{3}{4})$	0.9924	0.0076
  6	$[-\frac{1}{6},\frac{1}{3}]$	$(-\frac{1}{4},\frac{5}{12})$	0.6992	0.3008	18	$[-\frac{5}{6},\frac{2}{3}]$	$(-1,\frac{3}{4})$	0.9934	0.0066
  7	$[-\frac{1}{4},\frac{1}{3}]$	$(-\frac{1}{3},\frac{5}{12})$	0.7465	0.2535	19	$[-\frac{5}{6},\frac{3}{4}]$	$(-1,\frac{5}{6})$	0.9940	0.0060
  8	$[-\frac{1}{3},\frac{1}{3}]$	$(-\frac{5}{12},\frac{5}{12})$	0.8024	0.1976	20	$[-\frac{5}{6},\frac{5}{6}]$	$(-1,1)$	0.9992	0.0008
  9	$[-\frac{1}{3},\frac{5}{12}]$	$(-\frac{5}{12},\frac{1}{2})$	0.8371	0.1629	21	$[-\frac{5}{6},1]$	$(-1,1]$	0.9997	0.0003
  10	$[-\frac{1}{3},\frac{1}{2}]$	$(-\frac{5}{12},\frac{7}{12})$	0.9069	0.0931	22	$[-1,1]$	$[-1,1]$	1.0000	0.0000
  11	$[-\frac{5}{12},\frac{1}{2}]$	$(-\frac{1}{2},\frac{7}{12})$	0.9223	0.0777

  i is the iteration, MinCI is the smallest possible CI for the given CP and is always closed, MaxCI is the largest possible CI for the given CP and is open whenever is possible, and p_CI is the CP. The probabilities are given with four significant figures, however, the calculations and the actual values must be done with machine-like precision. The CI boundaries are given as irreducible fractions, and also for convenience, the integer fraction is reduced.

  and

  Table A8. Confidence intervals in increasing coverage for RR (data from Table A6).

  i	MinCI	MaxCI	${\mathit{p}}_{\mathit{CI}}$	$1-{\mathit{p}}_{\mathit{CI}}$	i	MinCI	MaxCI	${\mathit{p}}_{\mathit{CI}}$	$1-{\mathit{p}}_{\mathit{CI}}$	i	MinCI	MaxCI	${\mathit{p}}_{\mathit{CI}}$	$1-{\mathit{p}}_{\mathit{CI}}$
  0	$[\frac{4}{3},\frac{4}{3}]$	$(\frac{10}{9},\frac{5}{3})$	0.1563	0.8437	6	$[1,\frac{10}{3}]$	$(\frac{8}{9},1)$	0.3692	0.6308	12	$[\frac{1}{2},\frac{1}{0}]$	$(\frac{4}{9},\frac{1}{0}]$	0.8543	0.1457
  1	$[\frac{4}{3},\frac{5}{3}]$	$(\frac{10}{9},1)$	0.1598	0.8402	7	$[\frac{8}{9},\frac{10}{3}]$	$(\frac{5}{6},1)$	0.3730	0.6270	13	$[\frac{4}{9},\frac{1}{0}]$	$(\frac{1}{3},\frac{1}{0}]$	0.8698	0.1302
  2	$[\frac{4}{3},1]$	$(\frac{10}{9},\frac{8}{3})$	0.2527	0.7473	8	$[\frac{8}{9},1]$	$(\frac{5}{6},\frac{1}{0})$	0.3736	0.6264	14	$[\frac{1}{3},\frac{1}{0}]$	$(\frac{2}{9},\frac{1}{0}]$	0.9266	0.0734
  3	$[\frac{4}{3},\frac{8}{3}]$	$(\frac{10}{9},\frac{10}{3})$	0.2874	0.7126	9	$[\frac{8}{9},\frac{1}{0}]$	$(\frac{5}{6},\frac{1}{0}]$	0.6623	0.3377	15	$[\frac{2}{9},\frac{1}{0}]$	$(\frac{1}{6},\frac{1}{0}]$	0.9390	0.0610
  4	$[\frac{4}{3},\frac{10}{3}]$	$(\frac{10}{9},1)$	0.2943	0.7057	10	$[\frac{5}{6},\frac{1}{0}]$	$(\frac{2}{3},\frac{1}{0}]$	0.6623	0.3377	16	$[\frac{1}{6},\frac{1}{0}]$	$(0,\frac{1}{0}]$	0.9400	0.0600
  5	$[\frac{10}{9},\frac{10}{3}]$	$(1,1)$	0.2951	0.7049	11	$[\frac{2}{3},\frac{1}{0}]$	$(\frac{1}{2},\frac{1}{0}]$	0.8535	0.1465	17	$[0,\frac{1}{0}]$	$[0,\frac{1}{0}]$	1.0000	0.0000

  ).
• In order to gain a better understanding, another case is added here: the list of the successive expanding of the CI along with its CP. This is essentially the previous case with no stop limit and it is a more convenient arrangement of the data from Equation (15). Reviewing the described algorithm (Equation (16)):

  $$f_L\leftarrow \sum _{i_1\leftarrow L}^{i_0-1}{f}_{i_1},f_U\leftarrow \sum _{i_2\leftarrow i_0+1}^U{f}_{i_2},C{I}_{L,U}\leftarrow [h_L,h_U],p_{L,U}\leftarrow f_L+f_{i_0}+f_U$$

  (16)
  The example case listed in Table A4 (excess risk) is continued with the result of the above given recipe in Table 3. Examining entry 14 in Table 3, any interval between $[-\frac{7}{12},\frac{1}{2}]$ and $(-\frac{2}{3},\frac{7}{12})$ will produce the closest to (and greater than) the 95% (probability) coverage (which is about 95.40%, see entry 14 in Table 3); similarly, if one wants a 2% risk of being in error (a 98% CI), one should look to the entry 16 in Table 3; if 1% is the accepted level of risk of being in error, then the answer is on line 18 in Table 3 and, finally, for 1‰the answer is in line 21 in Table 3.
  Moving on to the example cases listed in the appendix in Table A5 (OR) and in Table A6 (RR), the results of the above given recipe are listed in

  Table A7. Confidence intervals in increasing coverage for OR (data from Table A5).

  i	MinCI	MaxCI	${\mathit{p}}_{\mathit{CI}}$	$1-{\mathit{p}}_{\mathit{CI}}$	i	MinCI	MaxCI	${\mathit{p}}_{\mathit{CI}}$	$1-{\mathit{p}}_{\mathit{CI}}$	i	MinCI	MaxCI	${\mathit{p}}_{\mathit{CI}}$	$1-{\mathit{p}}_{\mathit{CI}}$
  0	$[\frac{3}{2},\frac{3}{2}]$	$(1,\frac{5}{3})$	0.1389	0.8611	6	$[\frac{1}{5},\frac{3}{2}]$	$(\frac{1}{6},\frac{5}{3})$	0.4632	0.5368	12	$[0,1]$	$[0,1)$	0.6653	0.3347
  1	$[1,\frac{3}{2}]$	$(\frac{2}{3},\frac{5}{3})$	0.2130	0.7870	7	$[\frac{1}{6},\frac{3}{2}]$	$(\frac{1}{15},\frac{5}{3})$	0.4787	0.5213	13	$[0,1]$	$[0,1)$	0.6688	0.3312
  2	$[\frac{2}{3},\frac{3}{2}]$	$(\frac{3}{5},\frac{5}{3})$	0.2168	0.7832	8	$[\frac{1}{15},\frac{3}{2}]$	$(0,\frac{5}{3})$	0.4910	0.5090	14	$[0,1]$	$[0,1)$	0.7035	0.2965
  3	$[\frac{3}{5},\frac{3}{2}]$	$(\frac{1}{2},\frac{5}{3})$	0.3279	0.6721	9	$[0,\frac{3}{2}]$	$[0,\frac{5}{3})$	0.5546	0.4454	15	$[0,1]$	$[0,\frac{1}{0})$	0.7104	0.2896
  4	$[\frac{1}{2},\frac{3}{2}]$	$(\frac{1}{3},\frac{5}{3})$	0.3974	0.6026	10	$[0,\frac{5}{3}]$	$[0,1)$	0.5553	0.4447	16	$[0,\frac{1}{0}]$	$[0,\frac{1}{0}]$	1.0000	0.0000
  5	$[\frac{1}{3},\frac{3}{2}]$	$(\frac{1}{5},\frac{5}{3})$	0.4077	0.5923	11	$[0,1]$	$[0,1)$	0.5727	0.4273

  and Table A8, respectively (in the appendix too). For both (OR and RR of (x = 2, y = 1, m = 6, n = 4) configuration), the CI at no more than 5% risk of being in error is identical with the ones at 2%, 1% and 1‱ covering all possible cases, $[0,\infty ]$ (see entries 15 and 16 in Table A7 and 16 and 17 in Table A8).
• If someone wants to assess the probability of a no excess risk (ER, Equation (10), data arranged as in Equation (15)), that is the probability to pass through 0 trying to express the confidence; that is, in a much common language, the probability to be different from 0; then, the results given in the form of Equation (16) must be inspected. If ER ($\leftarrow \frac{xn-ym}{mn}$) is positive, then any positive $h_i$ ($\leftarrow \frac{un-vm}{mn}$) will satisfy the requirement and the associated probability is ${\sum }_{h_i>0}{f}_i$ and if ER < 0, then any $h_i$ < 0 will satisfy the requirement and the probability is ${\sum }_{h_i<0}{f}_i$. Inspecting Table 3, the last interval not containing 0 is at entry (iteration) 4, and its CP is 44.07%; the first interval containing 0 is at entry 5 in Table 3 and its CP is 51.47%; since the case of being 0 is the undesired, one should be excluded and, thus, the probability (for the excess risk, calculated at $\frac{1}{12}$, to be different from 0 for (x = 2, y = 1, m = 6, n = 4) configuration) is the one given at entry 4 in Table 3, 44.07%.
• Similarly, if someone wants to assess the probability of different odds (OR, Equation (11), data arranged as in Equation (15)), then again the results given in the form of Equation (16) must be inspected. If OR > 1, then any $h_i$ > 1 will satisfy the requirement and the associated probability is ${\sum }_{h_i>1}{f}_i$ and if OR < 1, then any $h_i$ < 1 will satisfy the requirement and the probability is ${\sum }_{h_i<1}{f}_i$. Inspecting Table A7, the last interval not containing 1 is at entry (iteration) 0, and its CP is 13.89%, thus giving a probability (for the OR, calculated at $\frac{3}{2}$, to be greater than 1 for (x = 2, y = 1, m = 6, n = 4) configuration) of 13.89%.
• The same is for a non-equal risk (RR, Equation (12), data arranged as in Equation (15)), the results given in Equation (16) providing the answer. If RR > 1, then any $h_i$ > 1 will satisfy the requirement and the probability is ${\sum }_{h_i>1}{f}_i$ and if RR < 1, then any $h_i$ < 1 will satisfy the requirement and the probability is ${\sum }_{h_i<1}{f}_i$. Inspecting Table A8, for the RR calculated at $\frac{4}{3}$, the last interval not containing 1 is at entry (iteration) 1, thus giving a probability (for RR > 1 for the (x = 2, y = 1, m = 6, n = 4) configuration) of 15.98%.

5.2. The Needed Algorithms for CIs of the Expressions of Two Proportions

Each step leading to the CIs for the expressions of two proportions is rigorously controlled by an algorithm. In this section, the five algorithms (Algorithms 1–5) are discussed and used to provide the results for x = 2, y = 1, m = 6 and n = 4 example. Thus,

• Algorithm 1 calculates the binomial probability mass function; the results of $\mathrm{B}(2,6,b1)$ and $\mathrm{B}(1,4,b2)$ (B$(z,o,r)$ in Algorithm 1) are in Table 2;
• Algorithm 2 calculates the ER (Equation (10), $\mathrm{E}\mathrm{R}(x,y,m,n,c,d)$), (OR (Equation (11), with $\mathrm{O}\mathrm{R}(x,y,m,n,c,d)$), and RR (Equation (12), with $\mathrm{R}\mathrm{R}(x,y,m,n,c,d)$) as irreducible fractions (with $\mathrm{F}\mathrm{R}(a,b,c,d)$) used in Algorithm 3 and later in succession with Algorithm 5; $\mathrm{E}\mathrm{R}(2,1,6,4,c,d)$ evaluates the ER with $\frac{c}{d}$ at $c=1$ and $d=12$; $\mathrm{O}\mathrm{R}(2,1,6,4,c,d)$ evaluates the OR with $\frac{c}{d}$ at $c=3$ and $d=2$; $\mathrm{R}\mathrm{R}(2,1,6,4,c,d)$ evaluates the RR with $\frac{c}{d}$ at $c=4$ and $d=3$;
• Algorithm 3 collects all possible drawings from the bivariate binomial distribution and pairs the probabilities with one (at the time) defined function operating on binomial proportion; the output of $\mathrm{B}2\mathrm{F}(x,y,m,n,\mathrm{F},p,g,h)$ is a $(m+1)\times (n+1)$ entries list (Equation (14)); the output of $\mathrm{B}2\mathrm{F}(2,1,6,4,\mathrm{F},p,g,h)$ for F ← ER is listed in Table A1, F ← OR in Table A2, and F ← RR in Table A3;
• Algorithm 4 with PMF_B2E(k,p,g,h) sorts and groups values and their associated probabilities; an output of it is formally given in Equation (15). Through it, data from Table A1 become the data in Table A4, from Table A2 to Table A5, and from Table A3 to Table A6;
• Algorithm 5 constructs CIs in increasing coverage and its output is formally given as Equation (16); the result of the evaluation of $\mathrm{C}\mathrm{I}\_\mathrm{B}2\mathrm{E}(nu,de,k,p,g,h)$ of the data in Table A4 (with $nu=1,de=12$) is listed in Table 3, of the data in Table A5 (with $nu=3,de=2$) in Table A7 and of the data in Table A6 (with $nu=4,de=3$) in Table A8.

In order to provide the desired output, the algorithms need to be chained together properly. Thus,

• So as to evaluate the ER (Equation (10)) for a certain ($x,y,m,n$) contingency (Figure 3), one should follow the sequence: $\mathrm{B}2\mathrm{F}(x,y,m,n,ER,p,g,h)$; $mn\leftarrow (m+1)·(n+1)$; $\mathrm{P}\mathrm{M}\mathrm{F}\_\mathrm{B}2\mathrm{E}(mn,p,g,h)$; $\mathrm{E}\mathrm{R}(x,y,m,n,r1,r2)$; $\mathrm{C}\mathrm{I}\_\mathrm{B}2\mathrm{E}(r1,r2,mn,p,g,h,ci)$ having then as output ($ci$) all the information required to construct Table 3;
• Similarly, in order to evaluate the OR (Equation (11)) for ($x,y,m,n$) contingency, one should follow the sequence: $\mathrm{B}2\mathrm{F}(x,y,m,n,OR,p,g,h)$; $mn\leftarrow (m+1)·(n+1)$; $\mathrm{P}\mathrm{M}\mathrm{F}\_\mathrm{B}2\mathrm{E}(mn,p,g,h)$; $\mathrm{O}\mathrm{R}(x,y,m,n,r1,r2)$; $\mathrm{C}\mathrm{I}\_\mathrm{B}2\mathrm{E}(r1,r2,mn,p,g,h,ci)$ having then as output ($ci$) all the information required to construct Table A7;
• Finally, in order to evaluate the RR (Equation (11)) for ($x,y,m,n$) the observed configuration, one should successively: $\mathrm{B}2\mathrm{F}(x,y,m,n,RR,p,g,h)$; $mn\leftarrow (m+1)·(n+1)$; $\mathrm{P}\mathrm{M}\mathrm{F}\_\mathrm{B}2\mathrm{E}(mn,p,g,h)$; $\mathrm{R}\mathrm{R}(x,y,m,n,r1,r2)$; $\mathrm{C}\mathrm{I}\_\mathrm{B}2\mathrm{E}(r1,r2,mn,p,g,h,ci)$, when the output ($ci$) is all the information required to construct Table A8.

5.3. Properties of the CIs

Some common beliefs about the CIs are that they should be larger in the middle (largest for values of x and y in the middle of $\{0,\dots ,m\}$ and $\{0,\dots ,n\}$, respectively), shorter at the ends (shortest for $x=0,m$ and $y=0,n$). In addition, they become shorter with the increase of the sample size (true if divided by the sample size), and their ends to be monotonic (slightly more complicated to keep the score here). As expected, most of the beliefs are true, but, even more so, some special properties were identified for the CIs. Thus, if $C{I}_{\mathrm{ER},\alpha }(x,y,m,n)$ = [$CI{L}_{\mathrm{ER},\alpha }(x,y,m,n)$, $CI{U}_{\mathrm{ER},\alpha }(x,y,m,n)$] (here the property is stated for the minimal CI, MinCI in Table 3 but is as well as true for the maximal one, MaxCI in Table 3), then $CI{L}_{\mathrm{ER},\alpha }(x,y,m,n)$ + $CI{U}_{\mathrm{ER},\alpha }(m-x,n-y,m,n)$ = 0 and $CI{U}_{\mathrm{ER},\alpha }(x,y,m,n)$ + $CI{L}_{\mathrm{ER},\alpha }(m-x,n-y,m,n)$ = 0 such that one could say that it is an odd function (in (m, n) modulo for (x, y)) and possessing in the middle some sort of antisymmetry too (see (3, 2) cell in Figure 5). This is an expected property since the excess risk itself possesses the same property ($\mathrm{ER}(x,y,m,n)$ + $\mathrm{ER}(m-x,n-y,m,n)$ = 0) and the construction of the CI (see the iteration for $p[i1-1]=p[i2+1]$ in Algorithm 5) gives the same chance to the bounds to be expanded.

A similar property possesses OR too (see Figure 6) but not RR (see Figure 7).

For the OR (see Figure 6) if [$\frac{CIL{0}_{\mathrm{OR},\alpha }(x,y,m,n)}{CIL{1}_{\mathrm{OR},\alpha }(x,y,m,n)}$, $\frac{CIU{0}_{\mathrm{OR},\alpha }(x,y,m,n)}{CIU{1}_{\mathrm{OR},\alpha }(x,y,m,n)}$] is the (minimal, see Table A7; the same applies for the maximal) CI for the OR in the ($x,y,m,n$) bivariate binomial experiment, then $CIL{0}_{\mathrm{OR},\alpha }(x,y,m,n)$ = $CIU{1}_{\mathrm{OR},\alpha }(m-x,m-y,m,n)$ and $CIL{1}_{\mathrm{OR},\alpha }(x,y,m,n)$ = $CIU{0}_{\mathrm{OR},\alpha }(m-x,m-y,m,n)$ as well as $CIU{0}_{\mathrm{OR},\alpha }(x,y,m,n)$ = $CIL{1}_{\mathrm{OR},\alpha }(m-x,n-y,m,n)$ and $CIU{1}_{\mathrm{OR},\alpha }(x,y,m,n)$ = $CIL{0}_{\mathrm{OR},\alpha }(m-x,n-y,m,n)$ such that the multiplication between $\frac{CIL{0}_{\mathrm{OR},\alpha }(x,y,m,n)}{CIL{1}_{\mathrm{OR},\alpha }(x,y,m,n)}$ and $\frac{CIU{0}_{\mathrm{OR},\alpha }(m-x,n-y,m,n)}{CIU{1}_{\mathrm{OR},\alpha }(m-x,n-y,m,n)}$ is formally 1 (and as well for $\frac{CIU{0}_{\mathrm{OR},\alpha }(x,y,m,n)}{CIU{1}_{\mathrm{OR},\alpha }(x,y,m,n)}$ and $\frac{CIL{0}_{\mathrm{OR},\alpha }(m-x,n-y,m,n)}{CIL{1}_{\mathrm{OR},\alpha }(m-x,n-y,m,n)}$). A simple check shows that this property is consistent with the property of OR: OR($x,y,m,n$) · OR($m-x,n-y,m,n$) = 1.

For small sample sizes (such as n = 4 and m = 6 from the selected example), the difference between the imposed level ($\alpha$) of being in error and the actual level of being in error can be big (see Figure 8, Figure 9 and Figure 10; but also cannot be bigger than the imposed level).

Figure 8, Figure 9 and Figure 10 show the real (actual) risks of being in error in constructing the CIs of ER (Figure 8), OR (Figure 9) and RR (Figure 10) for the selected example (CIs in Figure 5, Figure 6 and Figure 7). Since the process of constructing the CI is completely controlled (Algorithm 5), selecting the CI is easily set such that the actual error is the closest and smaller choice to the imposed level (see entry 14 in Table 3, entry 16 in Table A7 and entry 17 in Table A8). The difference between the imposed and actual level of being in error is dramatically diminished with the increase of the samples sizes. For comparison (with Figure 8), the next figure is given (Figure 11) in which it is visible that big departures (equal with the imposed level) occur only in the corners ($x=0,m$ and $y=0,n$) while the vast majority of the actual errors is in the [4.5%, 5.0%] range.

To be more precise, over 87% of the cases in the $x,y,m=60,n=50$ group of binomial experiment configurations depicted in Figure 11 have actual non-coverage probabilities in the [4.5%, 5.0%] range and about 95% of the cases have actual non-coverage probabilities in the [4.0%, 5.0%] range (2955 out of 3111 to be exact, the blue bullets in Figure 11).

There is no significant change in the shape of the actual coverage when changing from ER (Figure 11) to OR (Figure 12) and RR (Figure 13). For the ORs of $x,y,m=60,n=50$ group of binomial experiments depicted in Figure 12, the actual non-CPs fall in the [4.5%, 5.0%] range for 2775 cases (out of 3111, about 89%) and 2883 (over 92%) fall in the [4.0%, 5.0%] range (the blue bullets in Figure 12) while for the RRs (depicted in Figure 13), 2505 cases (about 80%) fall in the [4.5%, 5.0%] range and 2834 cases (about 91%) fall in the [4.0%, 5.0%] range (the blue bullets in Figure 13).

Considering that 220 cases from the $x,y,m=60,n=50$ group of binomial experiments have either x or y at the limits (0 or m and n, respectively) for which the non-coverage is very poor (only in 4 cases, CI at the limits for OR is in the [4.0%, 5.0%] range and in 62 for RR), for the rest of the cases ($x\in \{1,\dots ,m\}$, $y\in \{1,\dots ,n\}$), the accuracy is very good.

5.4. General Discussion

When dealing with finite populations, sampling strategy is foundational. Essentially, after extracting an individual (case) from a population, one can choose to put it back in the population (having a non-null likelihood to be extracted again in the future), or to throw it out. When the individual is returned to the population, a sampling with replacement is conducted, while when it is not returned to the population, a sampling without replacement is conducted.

Sampling with replacement is the one leading to the binomial distribution (by following a Bernoulli trial, counting successes vs. failures; given in Equation (1)). The alternative, sampling without replacement, leads to the hypergeometric distribution (Equation (17)).

$$f_H(k;n,K,N)\leftarrow \frac{\left(\genfrac{}{}{0pt}{}{K}{k}\right)\left(\genfrac{}{}{0pt}{}{N-K}{n-k}\right)}{\left(\genfrac{}{}{0pt}{}{N}{n}\right)}$$

(17)

When sampling without replacement (Equation (17)) the population (of size N in Equation (17)) is constantly diminished, so its size as well as the number of successes in the population (K) matters for the sample too (n—size; k—successes). Fisher has shown [68] that under sampling without replacement, the bivariate binomial experiment is transformed into a hypergeometric one; population is the union of the two samples. One alternative is to engage the sampling without replacement to extract from the population ($N\to m+n$) the first sample ($K\to m$) when sampling brings out from the total number of successes ($n\to x+y$) the successes from the (first) sample $k\to x$, when the probability (Fisher, $f_F$) to observe ($x,y,m,n$) experimental configuration is given in Equation (18).

$$f_F(x;y,m,n)\leftarrow f_H(x;x+y,m,m+n)=\frac{(x+y)!(m+n-x-y)!m!n!}{x!(m-x)!y!(n-y)!(m+n)!}$$

(18)

It must be noted that Equation (18) is symmetrical in $a=x,b=y,c=m-x,d=n-y$, which means that there exists an alternative way of engaging the sampling without replacement to extract from the population ($N\to m+n$) the same number (x) of successes. Sampling without replacement introduces a supplementary constraint to the experiment, such that dimensionality is reduced to 1 (Equation (19), Fisher Exact) and u ranges from 0 to $\min (m,n,x+y,m+n-x-y)$:

$$f_{FE}(u;x,y,m,n)\leftarrow \frac{(x+y)!(m+n-x-y)!m!n!}{u!(m-u)!(x+y-u)!(n+u-x-y)!(m+n)!}$$

(19)

One should notice that by substituting x with u and y with v in Equation (18), the right expression in Equation (18) becomes: $\frac{(u+v)!(m+n-u-v)!m!n!}{u!(m-u)!v!(n-v)!(m+n)!}$ and up to this point, two constrains have been used: the number of trials (draws) in the first experiment is m, and the number of trials (draws) in the second experiment is n; the supplementary constraint mentioned before is $x+y=u+v$, which when used leads to Equation (19). Fisher was the first to notice that the significance level is build up from the cases where the arrangement is as extreme as the observed arrangement, or more so. Some authors suggest that, when changing from evaluation of ER to OR or RR, to also change (see §1 in [61]) the study design and its associated statistical analysis (for Fisher Exact, see [69]) while others keep the bivariate binomial assumption even when some data are missing (see §2 in [70]). It is possible to argue that changing from m and n to $m+n$ does not change the nature of the phenomena (successions of Bernoulli trials) nor the independence between the samples; thus, Equation (9) is slightly changed replacing n and m with $m+n$ everywhere, but the increase of the cases leads to the increase in density of the points inside of the domain of the evaluated function and has as effect the increase of the CIs accuracy as well (for instance, for the case depicted in Figure 11, the increase is from 2995 to 3072 cases out of 3111 falling in the [4%, 5%] range).

Drawing with or without replacement from an infinite population should produce the same effect.

The proposed methodology (see Section 2, Section 5.1 and Section 5.2) is not significantly changed when the distribution is changed. Independent of the distribution, the methodology is:

• To collect all possible drawings and their associated probabilities (by using a possibly modified form of $\mathrm{B}2\mathrm{F}(x,y,m,n,\mathrm{F},p,g,h)$ from Algorithm 3) for one (at the time) defined function F (ER, Equation (10), $\mathrm{E}\mathrm{R}(x,y,m,n,c,d)$; OR, Equation (11), $\mathrm{O}\mathrm{R}(x,y,m,n,c,d)$; RR, Equation (12), $\mathrm{R}\mathrm{R}(x,y,m,n,c,d)$–all defined in Algorithm 2); the result is formally defined by Equation (14);
• Sort and group values and their associated probabilities (construct the probability mass function with PMF_B2E(k,p,g,h) from Algorithm 4); the result is formally defined by Equation (15);
• Construct the cumulative distribution function from which the CIs in increasing coverage (with $nu$ and $de$ from the evaluation of the selected (one of ER, OR, and RR) defined function (with $\mathrm{C}\mathrm{I}\_\mathrm{B}2\mathrm{E}(nu,de,k,p,g,h)$ from Algorithm 5); the result is formally defined by Equation (16).

In order to showcase the invariance of the strategy when the distribution is changed and to illustrate the alternating behavior of the actual non-CP, teeth-like plots, preferred by some authors (see Figure 4 in [8]), the actual non-CP calculated for the $m=9$ and $n=11$ group of experiment configurations (x from 0 to m and y from 0 to n) is depicted for three scenarios (Figure 14, Figure 15 and Figure 16) as function of the excess risk ($x/m-y/n$) when the imposed level of the non-CP was set to 5%.

Inspecting the Figure 14, Figure 15 and Figure 16, it is noticeable that the symmetry is kept in Figure 14 and Figure 16. Additionally, since the coverage and non-coverage probabilities are cumulative from different number of cases, the accuracy of the coverage is different and it is improved by the number of draws (about 4.2% for the data in Figure 15, $\mathrm{52,920}=10·12·{21}^2$ draws; about 3.66% for the data in Figure 14, 14,400 = ${10}^2·{12}^2$ draws; about 2.02% for the data in Figure 16, $1200=10·12·10$ draws).

The calculation of the exact CIs is of interest for other expressions of two proportions as well [71]. For a recent review of expressions involving two binomial proportions, one should consult [72], and when stratification is involved [73].

6. Conclusions

A strategy for providing exact confidence intervals for expressions of binomial proportions is elaborated (Figure 1). Algorithms for calculating the confidence intervals and/or p-values in the case of excess risk, odds ratio and relative risk are given (Algorithms 1–5) and exemplified. Extension of the findings for other distributions (continuous: Student t, Gaussian Normal; discrete: hypergeometric, joined samples binomial) is discussed.