Skip to content

fisher() returns a chi-square density on the wrong degrees of freedom, inverting the p-value #30

Description

@zhiaos

fisher() in methods/multilist.rs does not return a p-value. Two independent defects
combine so that the result is not merely inaccurate — it inverts, reporting the least
significant inputs as the most significant.

pub fn fisher(vals: &Vec<f64>) -> f64 {
    let k = vals.len();
    let pt = -2.0 * vals.iter().map(|x| x.ln()).sum::<f64>();
    let dist = statrs::distribution::ChiSquared::new(2_f64.powi(k as i32 - 1)).unwrap();
    dist.pdf(pt)
}

The statistic pt = -2 Σ ln(pᵢ) is correct. What follows is not.

1. .pdf() where an upper-tail probability is needed. Fisher's method wants
P(χ² ≥ pt), the survival function. A density is not a probability. For df > 2 the
chi-square density is zero at 0, rises to a mode, then decays — so it is small at both
ends. pt approaches 0 when every pᵢ approaches 1, which is the least significant
possible input, and that is exactly where the density is smallest.

2. 2^(k-1) degrees of freedom where 2k is needed. Exponential instead of linear:
k=3 gives 4 instead of 6, k=5 gives 16 instead of 10. Correct only at k=4, by coincidence.

Measured

input current correct
[0.999, 0.999, 0.999] 0.0015 1.0
[0.9, 0.9, 0.9] 0.115 0.996
[0.5, 0.5, 0.5] 0.130 0.655
[0.05, 0.05, 0.05] 0.00056 0.0063
[0.999] × 5 7.7e-21 1.0

Three p-values of 0.999 combine to 0.0015; five combine to 7.7e-21.

Reach

fisher() is pub, and MetaAnalysisMethod::Fisher is a public variant selected in
multilist_gsea (multilist.rs:96) and multilist_ora (multilist.rs:165). Any consumer of
webgestalt_lib choosing Fisher gets these values with no warning.

WebGestaltR is not currently affected. multiOraEnrichment.R:80 does pass "fisher",
so this runs on every multi-omics ORA analysis, but the R layer then discards the result and
recomputes the meta-p with poolr::stouffer (multiOraEnrichment.R:131, and
multiswGsea.R:178 for GSEA). The wrong value is computed and thrown away. That makes this
latent rather than live — and a hazard, since anyone removing the apparently redundant R
computation in favour of the Rust one would ship inverted p-values silently.

Suggested fix

pub fn fisher(vals: &Vec<f64>) -> f64 {
    let k = vals.len();
    let pt = -2.0 * vals.iter().map(|x| x.ln()).sum::<f64>();
    let dist = statrs::distribution::ChiSquared::new(2.0 * k as f64).unwrap();
    dist.sf(pt)          // or 1.0 - dist.cdf(pt) if `sf` is unavailable in this statrs version
}

sf is preferable to 1 - cdf — the latter loses precision in the far tail, which is where
combined p-values live.

Worth adding a doctest alongside the existing stouffer ones: combining [0.5, 0.5] should
give 0.5966 (χ² = 2.7726 on 4 df), and combining values near 1 should give a result near 1,
which is the property currently violated.

Metadata

Metadata

Assignees

Labels

bugSomething isn't working

Type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions