In the Background

Shortcuts in Fisher Exact P-values, part 2

Situation: We have a base assignment vector z\mathbf{z} of 1s and 0s, and we’re taking kk samples from the ntotal!n_\text{total}! permutations of z\mathbf{z}.

Symbols

  • ntreatn_\text{treat} : number of treated cases
  • nctrln_\text{ctrl} : number of control cases
  • ntotal=ntreat+nctrln_\text{total} = n_\text{treat} + n_\text{ctrl}
  • kk : number of draws of assignment vectors
  • Ncomb=(ntotalntreat)N_\text{comb} = \binom{n_\text{total}}{n_\text{treat}} : number of unique assignment vectors

Probability of any duplicates

Derivation

Goal: Probability that there are one or more duplicate assignment vectors in the kk samples. I.e., the probability that we can find at least 2 of the kk draws have the same pattern of 1s and 0s.

We’ll start with the probability that all draws are distinct assignment vectors, then invert it.

The first draw produces as assignment vector z1\mathbf{z}_1. The second permutation is drawn uniformly from the ntotal!n_\text{total}! permutations, and z2\mathbf{z}_2 (the assignment vector it produces) is drawn uniformly from the NcombN_\text{comb} assignment vectors. The probability that z2≠z1\mathbf{z}_2 \not = \mathbf{z}_1 is

Pr⁡(z2≠z1)=Ncomb−1Ncomb.\Pr(\mathbf{z}_2 \not = \mathbf{z}_1) = \frac{N_\text{comb} - 1}{N_\text{comb}}.

Likewise,

Pr⁡(z3∉{z1,z2}∣z2≠z1))=Ncomb−2Ncomb,\Pr(\mathbf{z}_3 \not \in \{\mathbf{z}_1, \mathbf{z}_2\} \mid \mathbf{z}_2 \not = \mathbf{z}_1)) = \frac{N_\text{comb} - 2}{N_\text{comb}},

and so on.

Therefore

P(zi all distinct )=1×Ncomb−1Ncomb×Ncomb−2Ncomb×⋯×Ncomb−(k−1)Ncomb=∏i=0k−1Ncomb−iNcomb\begin{aligned} P(\mathbf{z}_i \text{ all distinct }) &= 1 \times \frac{N_\text{comb} - 1}{N_\text{comb}} \times \frac{N_\text{comb} - 2}{N_\text{comb}} \times \cdots \times \frac{N_\text{comb} - (k-1)}{N_\text{comb}} \\ &= \prod_{i=0}^{k-1} \frac{N_\text{comb} - i}{N_\text{comb}} \end{aligned}

and

P(zi not all distinct )=1−∏i=0k−1(Ncomb−iNcomb)P(\mathbf{z}_i \text{ not all distinct }) = 1 - \prod_{i=0}^{k-1} \left(\frac{N_\text{comb} - i}{N_\text{comb}}\right)

Note: This is also the “birthday problem”.

Python code

from math import comb, log, exp  # comb requires version >= 3.8

def collision_prob(num_cases, num_treat, num_draws):
    num_combos = comb(num_cases, num_treat)

    if num_draws > num_combos:
        return 1.0

    # Compute the probability of no collisions
    # Using log and exp might help with large `num_combos`
    log_prob = sum(
        log(num_combos - j) - log(num_combos)
        for j in range(num_draws)
    )

    # Convert to probability of 1+ collisions
    return 1 - exp(log_prob)

Expected number of duplicates

Derivation

For each unique pattern i, the probability it appears at least once in k samples is

P(pattern i appears)=1−P(pattern i never appears)=1−(Ncomb−1Ncomb)k\begin{aligned} P(\text{pattern } i \text{ appears}) &= 1 - P(\text{pattern } i \text{ never appears}) \\ &= 1 - \left(\frac{N_\text{comb}-1}{N_\text{comb}}\right)^k \end{aligned}

By linearity of expectation,

E[# distinct]=∑i=1NcombP(pattern i appears)=Ncomb(1−(1−1Ncomb)k)\begin{aligned} E[\text{\# distinct}] &= \sum_{i=1}^{N_\text{comb}} P(\text{pattern } i \text{ appears}) \\ &= N_\text{comb} \left(1 - \left(1 - \frac{1}{N_\text{comb}}\right)^k\right) \end{aligned}

Then, since

# distinct+# duplicates=k,\text{\# distinct} + \text{\# duplicates} = k,

we have

E[# duplicates]=k−E[# distinct]=k−Ncomb(1−(1−1Ncomb)k)=k+Ncomb((1−1Ncomb)k−1)\begin{aligned} E[\text{\# duplicates}] &= k - E[\text{\# distinct}] \\ &= k - N_\text{comb} \left(1 - \left(1 - \frac{1}{N_\text{comb}}\right)^k\right) \\ &= k + N_\text{comb} \left(\left(1 - \frac{1}{N_\text{comb}}\right)^k - 1\right) \\ \end{aligned}

Expected fraction of duplicates

Derivation

Seems easy:

E[duplicates/k]=1+Ncombk((1−1Ncomb)k−1).E[\text{duplicates} / k] = 1 + \frac{N_\text{comb}}{k} \left(\left(1 - \frac{1}{N_\text{comb}}\right)^k - 1\right).

Unfortunately, this is difficult to compute when NcombN_\text{comb} is large.

A couple functions will help with accuracy

  • log1p(x) = log(1 + x)
  • expm1(x) = exp(x) - 1

Formula 1

Using log1p and expm1,

E[duplicates/k]=1+Ncombk((1−1Ncomb)k−1)=1+Ncombk(exp(log(1−1Ncomb)k)−1)=1+Ncombkexpm1(log(1−1Ncomb)k)=1+Ncombkexpm1(k⋅log(1−1Ncomb))=1+Ncombkexpm1(k⋅log1p(−1Ncomb))\begin{aligned} E[\text{duplicates} / k] &= 1 + \frac{N_\text{comb}}{k} \left(\left(1 - \frac{1}{N_\text{comb}}\right)^k - 1\right) \\ &= 1 + \frac{N_\text{comb}}{k} \left( \text{exp} \left( \text{log}\left(1 - \frac{1}{N_\text{comb}}\right)^k\right) - 1\right) \\ &= 1 + \frac{N_\text{comb}}{k} \text{expm1}\left(\text{log}\left(1 - \frac{1}{N_\text{comb}}\right)^k\right) \\ &= 1 + \frac{N_\text{comb}}{k} \text{expm1}\left(k \cdot \text{log}\left(1 - \frac{1}{N_\text{comb}}\right)\right) \\ &= 1 + \frac{N_\text{comb}}{k} \text{expm1}\left(k \cdot \text{log1p}\left(- \frac{1}{N_\text{comb}}\right)\right) \end{aligned}

Formula 2

For very large NcombN_\text{comb}, we’ll want to work with α=k/Ncomb\alpha = k / N_\text{comb}, and we’ll use

(1−1/Ncomb)k≈exp(−k/Ncomb)=exp(−α).(1 - 1 / N_\text{comb})^k \approx \text{exp}(-k / N_\text{comb}) = \text{exp}(- \alpha).

Then

E[duplicates/k]=1+Ncombk((1−1Ncomb)k−1)≈1+Ncombk(exp(−k/Ncomb)−1)=1+1α(exp(−α)−1)=1+expm1(−α)/α\begin{aligned} E[\text{duplicates} / k] &= 1 + \frac{N_\text{comb}}{k} \left(\left(1 - \frac{1}{N_\text{comb}}\right)^k - 1\right) \\ &\approx 1 + \frac{N_\text{comb}}{k} \left(\text{exp}(-k / N_\text{comb}) - 1\right) \\ &= 1 + \frac{1}{\alpha} \left(\text{exp}(-\alpha) - 1\right) \\ &= 1 + \text{expm1}(-\alpha) / \alpha \\ \end{aligned}

When NcombN_\text{comb} is very large, we’ll also want to use log(Ncomb)\text{log}(N_\text{comb}) in some intermediate calculations, via

Ncomb=exp(log(Ncomb)),N_\text{comb} = \text{exp}(\text{log}(N_\text{comb})),

and helper functions

  • gamma(x) = (x - 1)! when x is a positive integer, i.e., gamma(x + 1) = factorial(x)
  • gammaln(x) = log(gamma(x)).

Then

log(Ncomb)=log(ntotal!ntreat!⋅nctrl!)=gammaln(ntotal+1)−gammaln(ntreat+1)−gammaln(nctrl+1),\begin{aligned} \text{log}(N_\text{comb}) &= \text{log}\left( \frac{n_\text{total}!}{n_\text{treat}! \cdot n_\text{ctrl}!} \right) \\ &= \text{gammaln}(n_\text{total} + 1) - \text{gammaln}(n_\text{treat} + 1) - \text{gammaln}(n_\text{ctrl} + 1) \end{aligned},

and

α=exp(log(k)−log(Ncomb)).\alpha = \text{exp}\left(\text{log}(k) - \text{log}(N_\text{comb})\right).

Python code

from math import comb, log, exp, log1p, expm1
# comb requires version >= 3.8

from scipy.special import gammaln


def expected_duplicate_fraction(num_cases, num_treat, num_iter):
    """
    Expected fraction of draws that are duplicates: E[dup]/k

    Uses an approximation for large `num_cases`
    
    Parameters
    ----------
    num_cases : int - number of cases
    num_treat : int - number of 1's
    num_iter : int - number of draws
    
    Returns
    -------
    float - fraction of duplicates (between 0 and 1)
    """
    log_N = (
        gammaln(num_cases + 1) - 
        gammaln(num_treat + 1) - 
        gammaln(num_cases - num_treat + 1)
    )
    
    # Case 1: N is manageable
    if log_N < 700:
        N = comb(num_cases, num_treat)
        expected_distinct = -N * expm1(num_iter * log1p(-1/N))
        return 1 - expected_distinct / num_iter

    # Case 2: N is very large
    log_k = log(num_iter)
    alpha = exp(log_k - log_N)  # k/N

    if alpha < 1e-10:
        # For very small α, use Taylor series to avoid cancellation
        frac_dup = alpha / 2 - alpha**2 / 6 + alpha**3 / 24
    else:
        frac_dup = 1 + expm1(-alpha) / alpha
    
    return frac_dup