Overview of algorithm

A step by step overview of the algorithm and the helper functions it relies on

This article traces how the modified Fisher exact test is built in R, from the test frame through to the agreeing p-value and confidence interval. It assumes the statistical background, which is in the Background and comparison article; the derivations and proofs are in van der Meulen, Raymond and van der Meulen (2021).

The notation: \(u\) successes from \(m\) trials in Group A, \(v\) from \(n\) in Group B, with total \(T = u + v\).

	Group A: Success	Group A: Failure	Total
Group B: Success	..	..	\(v\)
Group B: Failure	..	..	\(v'\)
Total	\(u\)	\(u'\)	\(t\)

The basic premise

The modified test makes a single change to the randomised (UMPU) Fisher test: instead of always rejecting at a critical value (randomised) or never rejecting there (conservative), it uses one global threshold \(\gamma_0\) and rejects at a boundary only if that boundary’s randomisation probability \(\gamma\) exceeds \(\gamma_0\). The randomisation is taken away, and most of the boundary rejection mass is recovered.

The implementation chooses \(\gamma_0\) for a given null odds ratio, then inverts the resulting test for the p-value and confidence interval. The entry point is modified_fisher_exact_test():

modified_fisher_exact_test(u = 5, m = 12, v = 7, n = 11, odds_ratio = 1)

The steps below show what the above function does internally. Helper functions not exported to the namespace are prefixed by ., these can still be accessed through modifiedfisher::: (e.g. for troubleshooting purposes, or to visualise the steps with your own example as in this article).

Step 1: Build the test frame of critical values and randomisations

For a given null odds ratio \(\theta_0\) and each total \(T = 0, \ldots, m+n\), find the critical values \(c_1, c_2\) and admissible randomisation probabilities \(\gamma_1, \gamma_2\). This is construct_test_frame(), returning one row per \(T\) with columns c1, c2, gamma1, gamma2 (and the starting quantiles d1, d2).

With the critical values fixed, the two requirements (conditional size \(\alpha\) and unbiasedness) form a \(2\times2\) linear system solved for \((\gamma_1, \gamma_2)\) by .find_gamma12(). construct_test_frame() searches for the critical values giving an admissible solution (both \(\gamma\)’s in \([0, 1]\)), starting at the \(\alpha/2\) quantiles of the non-central hypergeometric distribution and spiralling outward in larger squares until one is found for every \(T\). At \(T = 0\) and \(T = m+n\) there is no room to randomise, so \(\gamma_1 = \gamma_2 = \alpha/2\). The frame holds \(2(m+n+1)\) randomisation probabilities in total, which matters in Step 4.

# Example 1 from the paper (m = 12, n = 11):
construct_test_frame(.odds_ratio = 1, .m = 12, .n = 11, .alpha = 0.05, .precision = 1e-3)

The first ten rows of this dataframe are given by:

#>    t c1 d1    gamma1 c2 d2    gamma2
#> 1  0  0  0 0.0250000  0  1 0.0250000
#> 2  1  0  0 0.0500000  1  2 0.0500000
#> 3  2  0  0 0.1100000  2  3 0.1000000
#> 4  3  0  0 0.2566667  3  4 0.2100000
#> 5  4  0  1 0.6416667  4  4 0.4666667
#> 6  5  1  1 0.1081481  4  5 0.0000337
#> 7  6  1  2 0.3630208  5  6 0.1892519
#> 8  7  1  2 0.9953571  6  6 0.5526948
#> 9  8  2  2 0.2761364  6  7 0.0543831
#> 10 9  2  3 0.8060101  7  8 0.3582323

Step 2: Define the rejection rule for a candidate \(\gamma_0\)

Given the frame and a candidate \(\gamma_0\), the verdict on any table \((u, T)\) is deterministic:

• \(u\) outside \([c_1, c_2]\): reject;
• \(u\) strictly inside \((c_1, c_2)\): accept;
• \(u = c_1\): reject iff \(\gamma_1 > \gamma_0\);
• \(u = c_2\): reject iff \(\gamma_2 > \gamma_0\);
• \(c_1 = c_2\) (a single boundary point): reject iff \(\gamma_1 + \gamma_2 > \gamma_0\).

This is .modified_reject(). Since the verdict depends only on the frame and \(\gamma_0\), the rejection region is built once into an \((m+1)\times(n+1)\) 0/1 matrix (.build_rejection_matrix()) and reused by local_size_modified(), power_modified(), and Step 4’s \(\gamma_0\) search. The \(\gamma_0\)-independent part (critical values and boundary masks) is cached separately (.build_rejection_base()), so Step 4’s search only repeats the cheap threshold comparison, not the full matrix construction.

Step 3: Measure the (worst case) size over the nuisance parameter \(\pi_1\)

A \(\gamma_0\) is acceptable only if the test’s true (unconditional) size stays below \(\alpha\). Conditional on \(T\) the test is size-\(\alpha\) by construction, but unconditionally the rejection probability depends on the unknown success rate \(\pi_1\). The real size is therefore a function of \(\pi_1\), and we need to guard against the ‘worst case’.

For a fixed \(\pi_1\), the rejection probability is a sum over all tables weighted by their two binomial probabilities, computed as the bilinear form \(p_u^\top R\, p_v\) with \(R\) the rejection matrix from Step 2 (local_size_modified()). The size is the maximum over \(\pi_1 \in (0, 1)\) (size_modified()).

Maximisation uses either "zoom" (grid, then re-grid finer around the peak, repeat; the default, and also the SAS macro’s approach) or "trust" (a trust-region optimiser with the analytic gradient .local_size_gradient_modified()). The size curve can be multi-peaked, so the safest check is to plot it via local_size_data = TRUE.

Both the zoom grid and the gradient are evaluated as batched matrix operations against R rather than as loops over individual grid points or tables, which is one of the reasons the R implementation outperforms the SAS macro (see the Reproducing the paper’s figures article).

testframe <- construct_test_frame(.odds_ratio = 1, .m = 12, .n = 11, .alpha = 0.05, .precision = 1e-3)

# Worst-case size at a candidate gamma0 of 0.05:
size_modified(.c = 0.05, .odds_ratio = 1, .m = 12, .n = 11, .df = testframe,
              .alpha = 0.05, .precision = 1e-3, .method = "zoom",
              .maze = 10, .zoom_iter = 6)

Step 4: Find the optimal \(\gamma_0\) by bisection

We still want the largest worst-case size to be less than \(\alpha\). The size only changes when \(\gamma_0\) crosses one of the \(2(m+n+1)\) actual \(\gamma\) values from Step 1, so rather than search a continuum, we sort those values and search among them.

Raising \(\gamma_0\) can only remove boundaries, so size is monotone in the sorted index, which makes a bisection appropriate. This is done in optimise_gamma0(): it sorts the pooled gamma1/gamma2 vector and bisects to the largest threshold whose size does not exceed \(\alpha\). The optimum is not a single number but a half-open interval between two adjacent sorted \(\gamma\)’s and any value in it gives the same test.

Since the test frame and the \(\gamma_0\)-independent part of the rejection matrix don’t change across this bisection, optimise_gamma0() builds each once and reuses them at every candidate; it also accepts an already-built test frame from the caller (.df) rather than rebuilding one, which matters because modified_fisher_exact_test() calls it repeatedly during the confidence interval and p-value searches.

optimise_gamma0(.odds_ratio = 1, .m = 12, .n = 11, .alpha = 0.05,
                .precision = 1e-3, .method = "zoom", .maze = 10, .zoom_iter = 6)

Step 5: Invert the test for the p-value and confidence interval

Steps 1 to 4 decide, at some given null odds ratio \(\theta_0\), whether the observed table is rejected (and the .accept() helper function ties it together: build the frame, optimise \(\gamma_0\), apply the rule). Confidence intervals and p-values are obtained by inverting that decision:

• Confidence interval: all \(\theta_0\) not rejected at level \(\alpha\). Its edges are found by bisection, starting from Woolf’s asymptotic limits and moving inward to the first rejection, to the requested precision.
• P-value: the smallest \(\alpha\) at which the table is rejected for the given null, found by bisecting on \(\alpha\).

Both invert the same test, so they agree by construction.

The main modified_fisher_exact_test() function runs the following pipeline:

construct_test_frame() → .modified_reject() / .build_rejection_matrix() → local_size_modified() / size_modified() → optimise_gamma0() → .accept()

The function returns an htest with $p.value,$estimate, $conf.int, the optimal $gamma0, the test frame, and optionally the size-versus-\(\pi_1\) dataframe:

result <- modified_fisher_exact_test(u = 5, m = 12, v = 7, n = 11, odds_ratio = 1, local_size_data = TRUE)

result$estimate
result$p.value
result$conf.int
result$local.size.data

Numerical limits and edge cases

Degenerate tables. If \(u = 0\) or \(u = m\) (all failures or all successes in one group), the test still runs but the confidence interval may be one-sided or very wide. The same applies when \(v = 0\) or \(v = n\). These are not errors; they reflect genuine uncertainty given the data.

Extreme odds ratios. Testing a null odds ratio far from the data (e.g. odds_ratio = 100 when the observed OR is near 1) is valid but may push the bisection searches to their limits. Set precision smaller (e.g. 1e-4) if the returned confidence limits look imprecise, and check that the returned p-value and confidence interval are still consistent.

Very small tables. The algorithm requires that \(m + n \geq 2\). For very small tables (e.g. \(m = n = 2\)), the test frame may have few rows and the optimised \(\gamma_0\) may be 0 (no borderline rejections), making the modified test identical to the conservative test. Check result$gamma0 and result$support.data to confirm the threshold is non-trivial.

Very large tables. Computation time grows with \((m+1)(n+1)\) because the rejection matrix has that many entries. The power = FALSE and conf_int = FALSE flags skip optional steps and reduce runtime when only the p-value is needed.

Grid density and size maximisation. The dominant cost in most calls is the maximisation of the local size over \(\pi_1\). This is controlled by the maze and zoom_iter arguments (for method = "zoom"). A denser \(\pi_1\) grid (equivalent to a larger \(N\) in the paper’s notation) gives a more reliable size estimate but takes proportionally longer: at \(m = 65\), \(n = 71\) the function takes 1-2 seconds at \(N = 20\) grid points and around 24 seconds at \(N = 300\). The default settings (maze = 10, zoom_iter = 6) are adequate for most purposes.

Precision. The default precision = 1e-3 is adequate for most purposes. Results are reported to three decimal places, matching the paper. For simulation studies or comparisons requiring finer agreement, use precision = 1e-4 or smaller, at the cost of additional computation.

Reference

van der Meulen EA, Raymond K, van der Meulen PJ (2021). Consistent Confidence Limits, P Values, and Power of the Non-Conservative, Size-α Modified Fisher Exact Test. Journal of Biostatistics and Biometric Applications 6(1):102.