cluster_composition's raw p-values stay anti-conservative after the dispersion correction

[timtreis]

The quasi-likelihood correction added in #88 fixes most of #82, but the raw p-values are still too small, and by an amount that depends on how many control wells there are.

Two reasons, both in _composition_test:

• The dispersion is estimated from the same control wells that define the composition every well is tested against (src/mantispy/tl/_heterogeneity.py:129 pools those controls into share, :170 takes their statistics as control_statistics). A control well contributes to the expected composition it is then compared with, so its chi-square is shrunk and the dispersion estimate is biased low.
• src/mantispy/tl/_heterogeneity.py:183 refers statistic / dispersion to a chi-square. The dispersion is itself estimated, which an F reference would account for and a chi-square does not.

Reproduction: draw every well of a plate from one composition with mild jitter, so nothing is a real hit, and count the wells called.

import warnings
import anndata as ad, numpy as np, pandas as pd
import mantispy as mt
from mantispy._core.schema import stamp

def wells(layout, n_controls):
    rows = [
        {"Metadata_Plate": "P1", "Metadata_Well": w, "Metadata_Perturbation": "DMSO" if i < n_controls else "pert",
         "Metadata_Control": i < n_controls, "leiden": str(c)}
        for i, (w, clusters) in enumerate(layout.items()) for c, n in clusters.items() for _ in range(n)
    ]
    obs = pd.DataFrame(rows, index=[str(i) for i in range(len(rows))])
    a = ad.AnnData(X=np.zeros((len(rows), 4), dtype=np.float32), obs=obs,
                   var=pd.DataFrame(index=[f"Cells_AreaShape_f{i}" for i in range(4)]))
    stamp(a, resolution="cell")
    return a

for n_controls in (8, 16, 32, 64):
    called = total = 0
    for trial in range(200):
        rng = np.random.default_rng(trial)
        shares = rng.dirichlet(np.full(4, 40.0), size=n_controls * 2)
        layout = {f"W{i:03d}": dict(enumerate(np.bincount(rng.choice(4, size=300, p=s), minlength=4)))
                  for i, s in enumerate(shares)}
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            comp = mt.tl.cluster_composition(wells(layout, n_controls))
        p = comp.uns["mantispy"]["composition_test"]["pvalue"].to_numpy()
        treated = ~comp.obs["Metadata_Control"].to_numpy(dtype=bool)
        called += int((p[treated] < 0.05).sum()); total += int(treated.sum())
    print(f"{n_controls:3d} controls: {called / total:.3f} called at p<0.05")

  8 controls: 0.138 called at p<0.05
 16 controls: 0.092 called at p<0.05
 32 controls: 0.070 called at p<0.05
 64 controls: 0.059 called at p<0.05

Expected 0.05 at every control count. The BH q-values are conservative on this null, so a caller who filters on qvalue is not affected; a caller who reads pvalue is.

The docstring already reports this residual rate in prose. The simulation behind it is not in the test suite, so the number cannot be re-derived from the repo.