[parashardhapola]

Summary

One study's list of changed genes rarely replicates in the next, so a result of several studies should say how well they agree next to what they find together. This adds combined p-values of several studies (Fisher, Stouffer with weights, Tippett, Wilkinson, Pearson, Mudholkar and George, Simes), one-sided or two-sided with the directions of the effects; Benjamini and Heller's partial conjunction p-value that an effect is found in at least u studies; Storey's proportion of true nulls (pi0, and pi1) as qvalue's pi0est computes it; and, for the studies of a CubeCollection, a replication analysis with sign agreement, heterogeneity, pi1 between studies and the genes found in k studies, with two figures. Each is checked against scipy, metap, qvalue, metafor and mpmath.

Changes

• cytepool.stats.combined_p_values.combine_p_values(p_values, *, method, weights, directions, rank): combines every feature's p-values (studies x features, NaN where a study did not test a feature) into one p-value of the global null hypothesis; returns CombinedPValues.
• cytepool.stats.combined_p_values.partial_conjunction(p_values, at_least, *, method, directions): the p-value that at least u of the studies have an effect, by Simes's, Fisher's or Stouffer's combination of the n - u + 1 largest p-values; returns CombinedPValues.
• cytepool.stats.combined_p_values.CombinedPValues: p-value, statistic, one-sided p-values up and down and the direction (with directions), and the number of studies per feature.
• cytepool.stats.multiple_testing.proportion_of_true_nulls(p, *, lambdas, method, smooth_df, smooth_log): pi0 as qvalue's pi0est, at one lambda, smoothed, or by the bootstrap; returns TrueNullProportion (pi0, pi1, lambdas, pi0_lambda, pi0_smooth).
• cytepool.replication.differential_expression_replication(collection, results, group, contrast, *, alpha, method, at_least, conjunction_method, weights, meta_method, null_proportion): reads one cell type and contrast through align_study_results and returns a ReplicationResult with the tables genes, study_pairs, found_in and studies, settings and provenance.
• cytepool.plotting.plot_study_agreement: a heatmap of pi1, or sign agreement, of every ordered pair of studies.
• cytepool.plotting.plot_replicated_genes: bars of the genes that exactly k studies find, split by whether every study finds them in the same direction.
• brent_minimum (private) takes an optional relative precision; a private smoothing spline in cytepool.stats.smoothing reproduces R's smooth.spline(x, y, df = df).
• Documentation: sections "Combined p-values", "Replicated in at least u studies", "The proportion of true nulls and pi1" and "Replication across the studies of a collection" in How the statistics work: Meta-analysis, a note on the Multiple testing page, API pages combined_p_values and replication, three rows of the capability matrix, 14 references (DOIs checked against Crossref) cited on the Bibliography page, the Validation page, the README's status table, planned list and limits, and the Overview.

Design notes

• Directions. The methods are defined for one-sided p-values. With directions (the signs of the effects), each two-sided p becomes p / 2 for an effect above zero and 1 - p / 2 below, as metap's two2one; both directions are combined and the p-value is min(1, 2 min(p_up, p_down)), the two-sided p-value of Z for Stouffer's method, so that effects in opposite directions cancel. The one-sided value near 0 is kept exactly (Stouffer's z is the quantile of p / 2 with the effect's sign, and every log is read from whichever of q and 1 - q is near 0), so p-values near 1e-300 keep their digits.
• Partial conjunction follows Benjamini and Heller: the n - u + 1 largest p-values, Stouffer unweighted (study weights attached to order statistics would break the validity argument). A gene that fewer than u studies tested has NaN and is left out of the Benjamini-Hochberg adjustment, a filter that does not depend on the p-values.
• One combination function serves both: combine_p_values is the partial conjunction with u = 1. Wilkinson's maximum (r = k) tests the global null, not that every study has an effect; the docs say so.
• Replication defaults: Stouffer weighted by the square roots of the samples each study's fit used, signed by the log2 fold changes; Simes's partial conjunction with u = 2; REML heterogeneity (NaN for genes of one study rather than metafor's Q = 0); pi1 with qvalue's defaults. Held-out prediction of case status (AUC) is not here: it needs per-sample scores of held-out studies, which belong with SumRank.
• Found by probing qvalue 2.38.0 (black box): it counts p-values at least lambda (the paper counts those above); its default lambda is R's seq(0.05, 0.95, 0.05), whose last value is 0.95 itself, not 0.05 + 18 x 0.05 = 0.9500000000000001; pi0 is the smoother's value at the largest lambda, capped at 1; the bootstrap is in closed form, the binomial variance of the count plus the squared distance of pi0(lambda) from the 10% quantile (type 7) of the pi0(lambda), first minimum (60 of 60 random cases identical); it drops NA.
• Found by probing R's smooth.spline (black box, documented outputs only): the penalty matrix (its keep.stuff output) integrates products of B-spline second derivatives with the weight 0.333 for the product of slopes, where the exact integral has 1/3; the ratio r uses the 3rd to the (nk - 3)th diagonal entries; the spar search is Brent's method with tolerance 1e-4 and relative precision 2e-8 (its trace output shows the steps and tol1); and the fit returned is that of the last spar evaluated, not of the best spar it reports. With these, pi0 agrees with qvalue to about 1e-12.
• Found by probing metap 1.1: it drops p-values of 0, and of 1 for sumz and logitp, and refuses one study; two2one rounds 1 - p / 2 to a double. scipy's combine_pvalues (BSD) was read for its conventions: Pearson's statistic is reported as -X, the logit's before scaling; no code was adapted.

Deliberate differences:

• p-values of 0 and 1 are combined as limits (Fisher's p-value with a 0 is 0; Stouffer's with both 0 and 1 is NaN), and one study gives its own p-value (except the logit, an approximation also for one study), where metap drops or refuses them.
• proportion_of_true_nulls gives an estimate where qvalue stops: where no p-value reaches the largest lambda (common for pi1 between studies: qvalue stopped for all 12 pairs of the validation with its defaults) and where the estimate is at most 0 (it is 0). The log smoother refuses a count of zero, as qvalue does. With more than 49 lambdas every lambda is a knot, where R keeps some (2.6e-5 relative for 96 lambdas).
• With directions, a two-sided p-value below 2.2e-308 loses digits when halved, and 5e-324 halves to 0.

Validation

validation/combined_p_values_and_replication.py, 54 checks: 47 pass, 7 deliberate differences, none failed (67.4 s, 444 MiB peak). Reference: scipy 1.18.1, mpmath 1.3.0; R 4.4.3, metap 1.1, qvalue 2.38.0, jsonlite 2.0.0; metafor 5.0.1.

Check family	Reference	Largest difference	Tolerance	Outcome
Fisher, Pearson, Tippett, weighted Stouffer, Mudholkar and George; 2, 5 and 20 studies, 0 to 15% missing	combine_pvalues	7.1e-15	1e-10	15 pass
Fisher, Stouffer, Tippett, Wilkinson, maximum, logit, Simes, as given and with directions (features with p >= 1e-5, where two2one keeps the digits)	metap, p.adjust	1.2e-12	1e-10	6 pass
Partial conjunctions, u = 1 to 3, Simes, Fisher, Stouffer, with and without directions	definition in R, with metap and p.adjust	7.8e-13	1e-10	6 pass
13 features of p-values from 5e-324 to 1 - 1e-15, every method, with and without directions	mpmath, 400 digits	2.0e-14	1e-10	1 pass
pi0: defaults, smooth.df = 5, smooth.log.pi0, bootstrap, one lambda, three grids, on 5 sets of 60 to 20,000 p-values	pi0est	1.5e-12	1e-10	8 pass
pi0 with a missing p-value; options both refuse	pi0est	6.7e-14	1e-10	4 pass
Replications of 3 simulated cohorts (2 cell types, 3 settings): genes, heterogeneity (REML, DerSimonian-Laird), pairs, pi1, genes found in k studies	metap, qvalue, p.adjust, rma.uni	2.0e-12	1e-10 (1e-6 for REML), 0 mismatches	5 pass
Figures: pi1 heatmap and bars as drawn	pi0est, R	0	1e-10, 0 mismatches	2 pass
metap at p = 0 and for one study; pi0est where no p-value reaches the largest lambda, where its estimate is at most 0, for 20 pairs of studies, and for 96 lambdas	metap, pi0est		n/a	7 deliberate differences

REF=/home/parashar/.cache/cytepool-reference
RSCRIPT=$REF/r-meta/bin/Rscript METAFOR_RSCRIPT=$REF/r-metafor/bin/Rscript \
    PYTHONPATH=. $REF/python/bin/python validation/combined_p_values_and_replication.py

Tests

tests/test_combined_p_values_and_replication.py, 90 tests, 1.1 s warm (22 s cold, compiling). They compare with fixtures the validation script writes (tests/data/combined_p_values_and_replication/, 72 KB: scipy, metap, mpmath, qvalue, and R's replication of three small cohorts rebuilt from their counts), with closed forms derived by hand (two studies, identical studies at Phi-bar(2), opposite signs, p-values of 0 and 1, one study, a gene in no study, the count of p-values at least lambda, p-values near 1e-300), and with invariants: the order of studies and of features, every direction turned around, a direction as the one-sided p-value of its sign, weights scaled together, Simes's partial conjunction growing with u, and the replication's columns against the public array functions applied to align_study_results. Messages of invalid arguments are tested; the masked-input inventory and the private-names list have the new callables.

Checks: the full suite on top of master, 4,571 passed, 18 skipped, also with bounds checks and with the rounding changed (OPENBLAS_CORETYPE=Sandybridge, 4 threads, NUMBA_CPU_NAME=generic); ruff check and ruff format clean; mypy clean; make -C docs api and make -C docs html (Node.js 24) build; tests/test_documentation.py 74 passed.