# Foundations Statistics Python

> Run scipy.stats/statsmodels tests (t-test, Mann-Whitney, ANOVA, chi-square, Fisher's exact, Pearson/Spearman) on expression/count/genotype data. Use when comparing groups, correcting p-values (FDR), or computing power/sample size in Python.

- Skill: `pavel-kravchenko/foundations-statistics-python` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/foundations-statistics-python`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/foundations-statistics-python/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/foundations-statistics-python

---


# Statistics with Python

## When to Use

- Comparing two or more groups of gene expression, protein abundance, or other continuous measurements (t-test, Mann-Whitney U, ANOVA, Kruskal-Wallis).
- Testing categorical/count associations: Hardy-Weinberg genotype counts (chi-square), SNP-disease 2x2 tables (Fisher's exact).
- Modeling variant allele counts, somatic mutation counts, or RNA-seq read counts with binomial/Poisson/negative-binomial distributions.
- Correlating two continuous variables (mRNA vs protein) or fitting linear/multiple regression with statsmodels.
- Correcting p-values across thousands of genes (Bonferroni/BH-FDR) or computing required sample size/power before an experiment.

## Version Compatibility

scipy ≥1.11, statsmodels ≥0.14, numpy ≥1.24, pandas ≥2.0, Python ≥3.10.

## Prerequisites

`pip install scipy statsmodels numpy pandas matplotlib`. Assumes familiarity with probability distributions (see `foundations-probability`) and basic pandas/numpy.

## Pitfalls

- **`ttest_ind` defaults to Student's t-test** (`equal_var=True`), which assumes equal variances. Pass `equal_var=False` explicitly for Welch's t-test — the safer default when variances may differ, as is common comparing tumor vs normal expression.
- **`mannwhitneyu` requires `alternative`**: always pass `alternative='two-sided'` explicitly — relying on the default risks version-dependent behavior.
- **`multipletests` returns a 4-tuple**: `(reject, pvals_corrected, alphacSidak, alphacBonf)` — unpack all four or index explicitly, don't assume only two values.
- **`pearsonr`/`spearmanr` return `(statistic, pvalue)`**: unpack both — `r = stats.pearsonr(x, y)` silently gives you a tuple, not a float.
- **Coordinate systems**: BED is 0-based half-open, VCF/GFF are 1-based inclusive — mixing them in downstream stats causes off-by-one errors.
- **Multiple testing**: never read raw p-values from a genome-wide scan; always apply FDR (Benjamini-Hochberg) before drawing conclusions.

## Core Distributions for Bioinformatics Data

**Goal:** model gene expression (normal), variant allele counts (binomial), somatic mutation counts (Poisson), and overdispersed RNA-seq counts (negative binomial) using `scipy.stats`'s consistent `.pdf/.pmf`, `.cdf`, `.sf`, `.ppf`, `.rvs` interface.

**Approach:** build the distribution object once with `stats.<dist>(params)`, then call methods on it — no need to re-derive formulas per distribution.

```python
import numpy as np
from scipy import stats

def variant_read_pvalue(alt_reads, coverage, expected_vaf=0.5):
    """P-value for observing >= alt_reads under a binomial null model.

    Use to flag a read count as unlikely under the expected variant
    allele frequency (e.g. 0.5 for a heterozygous germline SNP).
    """
    null_dist = stats.binom(n=coverage, p=expected_vaf)
    return null_dist.sf(alt_reads - 1)  # P(X >= alt_reads)


def mutation_burden_pvalue(observed_muts, background_rate):
    """P-value for a gene's somatic mutation count under a Poisson background model.

    background_rate: expected mutations per gene (e.g. from gene length x
    genome-wide mutation rate). Flags candidate driver genes.
    """
    bg_dist = stats.poisson(mu=background_rate)
    return bg_dist.sf(observed_muts - 1)


def nb_from_mean_dispersion(mu, alpha):
    """Convert bioinformatics (mean, dispersion) to scipy's nbinom(n, p).

    var = mu + alpha * mu**2 (DESeq2/edgeR convention).
    """
    n = 1 / alpha
    p = 1 / (1 + alpha * mu)
    return stats.nbinom(n=n, p=p)


# Example: is 6 alt reads more consistent with germline (p=0.5) or somatic (p=0.2)?
germ_dist = stats.binom(n=30, p=0.5)
somatic_dist = stats.binom(n=30, p=0.2)
likelihood_ratio = somatic_dist.pmf(6) / germ_dist.pmf(6)
print(f"P(alt=6) is {likelihood_ratio:.1f}x more likely under the somatic model")
```

## Choosing and Running the Right Hypothesis Test

**Goal:** pick the correct test for the study design (parametric vs non-parametric, paired vs independent, two groups vs many) and run it correctly.

**Approach:** check assumptions first (`stats.shapiro` for normality, `stats.levene` for equal variance), then dispatch to the matching scipy/statsmodels function.

```python
import numpy as np
from scipy import stats
from statsmodels.stats.multicomp import pairwise_tukeyhsd

def check_assumptions(group_a, group_b, alpha=0.05):
    """Test normality (Shapiro-Wilk) and equal variance (Levene) for two groups.

    Returns a dict with booleans used to choose the downstream test.
    """
    _, p_norm_a = stats.shapiro(group_a)
    _, p_norm_b = stats.shapiro(group_b)
    _, p_var = stats.levene(group_a, group_b)
    return {
        "normal": (p_norm_a > alpha) and (p_norm_b > alpha),
        "equal_var": p_var > alpha,
    }


def compare_two_groups(group_a, group_b, alpha=0.05):
    """Auto-select t-test (Student/Welch) or Mann-Whitney U based on assumptions.

    Returns (test_name, statistic, p_value).
    """
    assumptions = check_assumptions(group_a, group_b, alpha)
    if assumptions["normal"]:
        stat, p = stats.ttest_ind(group_a, group_b, equal_var=assumptions["equal_var"])
        name = "Student's t-test" if assumptions["equal_var"] else "Welch's t-test"
    else:
        stat, p = stats.mannwhitneyu(group_a, group_b, alternative="two-sided")
        name = "Mann-Whitney U"
    return name, stat, p


# Multi-group comparison: ANOVA -> Tukey HSD post-hoc if significant
liver, kidney, brain = (np.random.normal(m, 1.2, 20) for m in (8.0, 8.3, 7.8))
f_stat, p_anova = stats.f_oneway(liver, kidney, brain)
if p_anova < 0.05:
    labels = ["liver"] * 20 + ["kidney"] * 20 + ["brain"] * 20
    tukey = pairwise_tukeyhsd(np.concatenate([liver, kidney, brain]), labels, alpha=0.05)
    print(tukey.summary())

# Categorical association: chi-square (large n) or Fisher's exact (small n / 2x2)
observed = np.array([[38, 12], [22, 28]])   # [disease, control] x [alt allele, no alt allele]
odds_ratio, p_fisher = stats.fisher_exact(observed, alternative="two-sided")
print(f"Odds ratio: {odds_ratio:.3f}, p={p_fisher:.4f}")
```

## Correlation, Regression, and Multiple Testing Correction

**Goal:** quantify a relationship between two continuous variables and correct p-values when testing many features (e.g. every gene) at once.

**Approach:** use `pearsonr`/`spearmanr` for a single pair, `statsmodels.formula.api.ols` for regression with multiple predictors, and `multipletests` for FDR correction across a p-value array.

```python
import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf
from statsmodels.stats.multitest import multipletests

def correlate_and_regress(mrna, protein):
    """Pearson + Spearman correlation, then simple OLS: protein ~ mrna.

    Spearman is more robust when protein data contain outliers.
    """
    r_pearson, p_pearson = stats.pearsonr(mrna, protein)
    r_spearman, p_spearman = stats.spearmanr(mrna, protein)
    df = pd.DataFrame({"mrna": mrna, "protein": protein})
    model = smf.ols("protein ~ mrna", data=df).fit()
    return {
        "pearson_r": r_pearson, "pearson_p": p_pearson,
        "spearman_r": r_spearman, "spearman_p": p_spearman,
        "ols_r_squared": model.rsquared, "ols_slope_p": model.pvalues["mrna"],
    }


def fdr_correct(pvalues, alpha=0.05, method="fdr_bh"):
    """Benjamini-Hochberg FDR correction across many p-values (e.g. per-gene DE tests).

    Returns a DataFrame with the original and adjusted p-values plus a
    significance flag — safe to sort/filter directly.
    """
    reject, pvals_corrected, _, _ = multipletests(pvalues, alpha=alpha, method=method)
    return pd.DataFrame({
        "pvalue": pvalues,
        "padj": pvals_corrected,
        "significant": reject,
    })


pvals = np.random.beta(0.5, 5, size=1000)  # simulated per-gene p-values
results = fdr_correct(pvals)
print(f"{results['significant'].sum()} genes significant at FDR 0.05")
```

## See Also

- `foundations-probability` — probability distributions and the theory behind the tests here.
- `foundations-biostatistics-fundamentals` — conceptual grounding for hypothesis testing.
- `bio-experimental-design-multiple-testing` — deeper coverage of FDR/Bonferroni strategies.
- `bio-experimental-design-power-analysis` — sample-size and power calculations with `statsmodels.stats.power`.