# Foundations Biostatistics Fundamentals

> Run t-test/Mann-Whitney/ANOVA/chi-square tests, BH-FDR correction, and power analysis in SciPy, statsmodels, R. Use when comparing groups, interpreting p-values/CIs, correcting many gene-level tests, or sizing an experiment.

- Skill: `pavel-kravchenko/foundations-biostatistics-fundamentals` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/foundations-biostatistics-fundamentals`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/foundations-biostatistics-fundamentals/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-biostatistics-fundamentals

---


# Biostatistics Fundamentals

## When to Use
- Deciding whether an observed difference (expression, protein level, allele frequency) between two or more groups is likely real or noise.
- Choosing between a parametric test (t-test, ANOVA, Pearson) and its non-parametric counterpart (Mann-Whitney, Kruskal-Wallis, Spearman).
- Correcting p-values across thousands of simultaneous tests (genes, variants, metabolites) with Bonferroni or Benjamini-Hochberg FDR.
- Computing a confidence interval around a mean/proportion, or estimating the sample size needed to detect an effect (power analysis).
- Sanity-checking a downstream tool's output (DESeq2, edgeR, GWAS) by understanding what the p-value and effect size actually mean.

## Version Compatibility
Python ≥3.10, NumPy ≥1.24, SciPy ≥1.11, statsmodels ≥0.14, matplotlib ≥3.7. R ≥4.3 with base `stats` (no extra packages needed for the tests shown here).

## Prerequisites
- `pip install numpy scipy statsmodels matplotlib`
- Comfortable with arrays/DataFrames; no prior stats background assumed, but you should know your experimental design (paired vs. independent groups, number of comparisons).

## Descriptive Statistics and Confidence Intervals

**Goal:** Summarize a sample's center/spread and quantify uncertainty in the estimated mean.
**Approach:** Report the median/IQR alongside the mean/SD (robust to outliers), and use the SEM with a t-distribution critical value for the confidence interval — not the normal distribution, unless n is large.

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


def descriptive_summary(x: np.ndarray, confidence: float = 0.95) -> dict:
    """Return center, spread, and a confidence interval for the mean of x.

    Uses the t-distribution (not z) since the population SD is unknown,
    which is almost always true for biological samples.
    """
    x = np.asarray(x, dtype=float)
    n = len(x)
    mean, sd = np.mean(x), np.std(x, ddof=1)
    sem = sd / np.sqrt(n)
    q25, q75 = np.percentile(x, [25, 75])
    t_crit = stats.t.ppf((1 + confidence) / 2, df=n - 1)
    ci = (mean - t_crit * sem, mean + t_crit * sem)
    return {
        "n": n, "mean": mean, "median": np.median(x), "sd": sd,
        "sem": sem, "iqr": q75 - q25, "ci": ci,
    }


expression = np.append(np.random.default_rng(42).normal(10, 2, 100), [20, 22, 25])  # outliers
print(descriptive_summary(expression))
```

## Choosing and Running a Hypothesis Test

**Goal:** Test whether two (or more) groups differ, using the test that matches the study design and data shape.
**Approach:** Check normality visually (Q-Q plot) plus Shapiro-Wilk as a rough guide (unreliable at very large/small n), match paired vs. independent design, then pick parametric or rank-based.

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


def compare_two_groups(a: np.ndarray, b: np.ndarray, paired: bool = False, alpha: float = 0.05) -> dict:
    """Compare two groups with the appropriate t-test or its non-parametric counterpart.

    Runs Shapiro-Wilk on both groups; if either rejects normality (p < alpha),
    falls back to Mann-Whitney U (independent) or Wilcoxon signed-rank (paired).
    """
    a, b = np.asarray(a, float), np.asarray(b, float)
    normal_a = stats.shapiro(a).pvalue > alpha
    normal_b = stats.shapiro(b).pvalue > alpha

    if normal_a and normal_b:
        if paired:
            stat, p = stats.ttest_rel(a, b)
            test = "paired t-test"
        else:
            stat, p = stats.ttest_ind(a, b)
            test = "independent t-test"
    else:
        if paired:
            stat, p = stats.wilcoxon(a, b)
            test = "Wilcoxon signed-rank"
        else:
            stat, p = stats.mannwhitneyu(a, b, alternative="two-sided")
            test = "Mann-Whitney U"

    pooled_sd = np.sqrt((np.var(a, ddof=1) + np.var(b, ddof=1)) / 2)
    cohens_d = (np.mean(a) - np.mean(b)) / pooled_sd
    return {"test": test, "statistic": stat, "p_value": p, "cohens_d": cohens_d}


rng = np.random.default_rng(42)
normal_tissue = rng.normal(10, 2, 30)
tumor_tissue = rng.normal(12, 2.5, 30)
print(compare_two_groups(normal_tissue, tumor_tissue))
```

For 3+ groups use `stats.f_oneway` (ANOVA, parametric) or `stats.kruskal` (Kruskal-Wallis, non-parametric). For categorical 2x2 tables use `stats.chi2_contingency` (large n) or `stats.fisher_exact` (small n, gives an odds ratio).

## Multiple Testing Correction and Power Analysis

**Goal:** Control the false discovery rate when testing thousands of genes/variants at once, and determine the sample size needed before running the experiment.
**Approach:** Never read raw p-values genome-wide — always correct with BH/FDR (less conservative than Bonferroni, standard for genomics). Before collecting data, use power analysis to justify replicate counts.

```python
import numpy as np
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.power import TTestIndPower


def correct_pvalues(pvalues: np.ndarray, alpha: float = 0.05, method: str = "fdr_bh") -> dict:
    """Adjust p-values for multiple comparisons and report how many pass."""
    reject, adj_p, _, _ = multipletests(pvalues, alpha=alpha, method=method)
    return {"n_significant": int(reject.sum()), "adjusted_pvalues": adj_p, "reject": reject}


def samples_needed(effect_size: float, power: float = 0.80, alpha: float = 0.05) -> int:
    """Minimum samples per group (two-sided independent t-test) for the given power."""
    n = TTestIndPower().solve_power(effect_size=effect_size, alpha=alpha, power=power, alternative="two-sided")
    return int(np.ceil(n))


rng = np.random.default_rng(42)
pvals = np.concatenate([rng.uniform(0, 1, 9800), rng.beta(0.3, 10, 200)])  # 200 truly DE genes
result = correct_pvalues(pvals)
print(f"Significant after BH-FDR: {result['n_significant']} / {len(pvals)}")
print(f"n per group for medium effect (d=0.5), 80% power: {samples_needed(0.5)}")
```

Equivalent workflow in R (base `stats`, no extra packages):

```r
## Two-group comparison with normality check, then p-value correction
group_a <- rnorm(30, mean = 10, sd = 2)
group_b <- rnorm(30, mean = 12, sd = 2.5)

shapiro.test(group_a)$p.value  ## if < 0.05, prefer wilcox.test below
shapiro.test(group_b)$p.value

t_result <- t.test(group_a, group_b)          ## Welch's t-test by default
wilcox_result <- wilcox.test(group_a, group_b)

cat("t-test p =", t_result$p.value, "\n")
cat("Mann-Whitney p =", wilcox_result$p.value, "\n")

## Multiple testing correction across many genes
pvals <- c(runif(9800), rbeta(200, 0.3, 10))
padj <- p.adjust(pvals, method = "BH")
cat("Significant after BH-FDR:", sum(padj < 0.05), "\n")

## Power analysis: samples per group for medium effect size (Cohen's d = 0.5)
power.t.test(delta = 0.5, sd = 1, power = 0.80, sig.level = 0.05)$n
```

## Pitfalls
- **P-value is not the probability that H0 is true**: it is P(data this extreme | H0 true). Confusing these two is the most common statistics error in papers.
- **Statistical significance ≠ biological significance**: with n = 10,000, a 0.1% fold-change will be "significant." Always report effect size (Cohen's d, log2FC) alongside p-values.
- **Multiple testing explosion**: testing 20,000 genes at alpha = 0.05 expects 1,000 false positives even with zero true effects — always apply BH/FDR correction in genomics.
- **Paired vs. unpaired mismatch**: running an independent t-test on before/after data from the same patients discards information and loses power; match the test to the design.
- **Shapiro-Wilk is not a decision oracle**: it over-rejects normality at large n and under-detects it at small n. Look at a Q-Q plot and use domain knowledge, not just the p-value.
- **R's `t.test` defaults to Welch's (unequal variance)**, while `scipy.stats.ttest_ind` defaults to equal-variance (Student's) unless you pass `equal_var=False` — these can give different p-values on the same data.

## See Also
- `bio-experimental-design-power-analysis` — deeper sample-size and power calculations for specific study designs.
- `bio-experimental-design-multiple-testing` — FDR/Bonferroni strategies beyond the basics shown here.
- `bio-differential-expression-de-results` — applying these tests to real RNA-seq count data (DESeq2/edgeR output).
- `bio-population-genetics-association-testing` — hypothesis testing applied to genotype-phenotype association (GWAS).