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.

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.

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.

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):

## 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).