permfdr

A Python-only toolkit for two recurring problems in permutation testing, implementing (porting) published methods — no invented statistics.

1. Tail approximation — smooth, strictly-positive p-values below the empirical floor 1/(Np+1) when permutations are expensive. Faithful port of permApprox (constrained generalized-Pareto tail model; Peschel et al. 2026).
2. FDR on discrete p-values — multiple-testing procedures that exploit the known discrete null of permutation p-values instead of assuming a continuous Uniform(0,1): DS-FDR (Jiang 2017), DBH (Döhler 2018), BH+ (Chen 2020), plus BH/BY/Storey/mid-p baselines.

See REVIEW.md for the method survey and GUIDE.md for a "when to use what" decision guide.

Install

Work inside the dedicated conda env (permfdr); this repo was built and tested against it (Python 3.11, OpenBLAS numpy/scipy):

python -m pip install -e .        # from the permfdr/ directory
pytest -q                         # 69 tests

Quick start

import numpy as np
from permfdr.permutation import permutation_test

# data: (n_features, n_samples), labels: (n_samples,)
obs, perm = permutation_test(data, labels, statistic="t", n_perm=2000, seed=0)
# obs: (n_features,)   perm: (n_perm, n_features)   -- permutations along axis 0

Problem 1 — sub-floor p-values (constrained GPD)

from permfdr.gpd import perm_approx

res = perm_approx(obs, perm, method="gpd", alternative="two_sided",
                  sample_size=n_smaller_group)   # SLLS margin needs a size proxy
res.p_values      # hybrid: empirical in the bulk, constrained-GPD in the tail
res.p_empirical   # raw empirical p (pinned at 1/(B+1) in the tail)
res.method_used   # 'empirical' or 'gpd' per test

perm_approx screens tests with p_empirical < 0.1, fits a GPD to the tail with a support constraint + adaptive safety margin, and guarantees p > 0.

Problem 2 — FDR across tests

from permfdr.fdr import (permutation_support, bh, dbh, bh_plus, dsfdr)

# bridge: observed p-values + their KNOWN discrete supports
pvals, supports = permutation_support(obs, perm, alternative="two_sided")

bh(pvals, alpha=0.05)["reject"]                       # reference
dbh(pvals, supports, alpha=0.05, direction="su")      # discrete BH (Döhler)
bh_plus(pvals, supports, alpha=0.05, adaptive=True)   # adaptive BH+ (Chen)
dsfdr(obs, perm, alpha=0.05)                           # permutation FDR (Jiang)

For count data use the exact discrete tests instead of the permutation bridge:

from permfdr.fdr import binomial_pvalues_support, fisher_pvalues_support, dbh
pvals, midp, supports = binomial_pvalues_support(counts, sizes)
dbh(pvals, supports, 0.05, "su")["n_reject"]

What is proven (tests + studies)

• GPD estimators recover known parameters and match scipy.genpareto; the constrained tail p-value is accurate (p_GPD/p_exact ≈ 1) and never zero, while the unconstrained (Knijnenburg) baseline produces zeros (studies/gpd_validation.py).
• DS-FDR matches the biocore/dsFDR reference exactly on identical inputs.
• Discrete p-values match scipy.binomtest / fisher_exact exactly; BH reproduces the R papers' published counts on Lister (420) and Gilbert (16) when fed R's own p-values (studies/reproduce_counts.py).
• FDR control + power: all methods control the FDR at q, with the power ordering DS-FDR ≥ DBH-SU ≥ aBH+ ≥ BH on discrete data (studies/fdr_control_sim.py).
• Real data, known signal: breast-cancer (permutation-p vs analytic-t correlation = 1.000; signal survives correction) and Golub leukemia (~1200 DE genes) — studies/breast_cancer_study.py, studies/golub_study.py.

Run a study: python studies/<name>.py (figures land in studies/figures/).

Layout

src/permfdr/
  permutation.py  empirical_p.py  normalize.py   # run tests, compute p, normalize
  gpd/            estimators, threshold, constrained (SLLS), approx, baselines, gof
  fdr/            support (bridge), bh, dsfdr, discrete_bh, chen, discrete_tests
  external/dsfdr_ref.py                           # vendored biocore core (BSD)
  datasets.py                                     # local + downloadable datasets
tests/   studies/   data/

License

GPL-3.0-or-later — the package contains Python translations of GPL-3 R code (permApprox, DiscreteFDR) with per-module attribution. The vendored permfdr.external.dsfdr_ref keeps its original BSD-3-Clause license.