Practical Biostatistics in 30 Days

From p-values to pipelines — a structured journey through the statistics every biologist actually needs.

You have the data. Thousands of gene expression measurements. Hundreds of patient outcomes. Millions of variants. You know the biology. You understand the experiment. But when it comes time to choose a statistical test, set a significance threshold, or interpret a confidence interval, the ground shifts under your feet.

This book fixes that. In 30 days, you will go from statistical anxiety to statistical fluency — not by memorizing formulas, but by solving real biological problems with real data. Every test you learn has a reason. Every formula has a story. Every p-value has a context.

And you will do it all in BioLang, a language with over 400 statistical builtins that lets you express an entire analysis — from data loading to hypothesis testing to publication-quality visualization — in a handful of readable, pipe-chained lines.

Who This Book Is For

This book is for anyone who works with biological data and needs to make sound statistical decisions. You might be:

• A biologist who dreads the statistics section. You can design elegant experiments, but when the reviewer asks why you used a t-test instead of a Mann-Whitney U, you panic. You have tried statistics textbooks, but they are full of coin flips and card games when you need differential expression and survival curves. This book teaches statistics through biology, using datasets and questions you actually care about.

• A developer entering biotech. You can write production code and build data pipelines, but you do not know the difference between a parametric and a non-parametric test. You have heard that bioinformatics requires “statistical thinking,” but nobody has explained what that means in practice. This book gives you the statistical intuition alongside the implementation, so you understand why you are computing a fold change, not just how.

• A graduate student facing qualifying exams. Your program expects fluency in biostatistics, but your coursework is a blur of Greek letters and proof sketches. You need a practical guide that connects the math to the biology and shows you how to actually run these tests on real data. This book builds that bridge in 30 structured days.

• A clinical researcher designing or analyzing studies. You work with patient cohorts, treatment outcomes, and survival data. You need to choose the right test, compute adequate sample sizes, and report results that satisfy both statisticians and regulatory reviewers. This book covers clinical biostatistics end to end — from power analysis through Cox proportional hazards.

No matter which category you fall into, you share one thing: you want statistical skills that solve real problems, not abstract exercises. Every day in this book produces an analysis you can adapt to your own data.

Your Path Through the Book

The first week builds foundations for every reader, but your starting point may differ. Here is which days to prioritize based on your background:

Complete beginner? That is completely fine. Day 1 starts with a single concept — what is a distribution? — and builds from there. No calculus. No linear algebra. No prior programming experience beyond basic BioLang familiarity. If you can type bl run script.bl, you are ready.

What You Will Learn

Over 30 days, you will go from statistical uncertainty to being able to:

• Describe and summarize any biological dataset (distributions, central tendency, spread, outliers)
• Choose the correct statistical test for any experimental design
• Perform and interpret t-tests, ANOVA, chi-square tests, and their non-parametric alternatives
• Run linear and logistic regression on biological data
• Analyze time-to-event data with Kaplan-Meier curves and Cox proportional hazards models
• Reduce high-dimensional data with PCA and interpret biplots
• Cluster samples and genes using hierarchical and k-means methods
• Correct for multiple testing with Bonferroni, Benjamini-Hochberg, and permutation approaches
• Compute effect sizes, confidence intervals, and statistical power
• Design experiments with proper sample size calculations
• Build volcano plots, Manhattan plots, Q-Q plots, and forest plots
• Apply Bayesian reasoning to biological problems
• Complete three capstone analyses that mirror real research publications

You will learn all of this in BioLang, which provides dedicated builtins for every test and method. But you will not be locked in. Every day includes comparison examples in Python (scipy/statsmodels) and R (base stats/survival/ggplot2), so you can translate your skills to any environment.

How This Book Is Structured

The book is organized into four weeks plus capstone projects:

Each day follows the same structure:

1. The Problem — a vivid scenario that shows why you need today’s method. A researcher staring at ambiguous results. A clinician choosing between treatments. A graduate student defending a finding.
2. What Is [Topic]? — a plain-language explanation of the statistical concept, free of jargon. If your collaborator asked “what is a p-value?” at a coffee shop, this is how you would explain it.
3. Core Concepts — the ideas, assumptions, and mechanics, presented with tables, diagrams, and worked examples. Formulas appear when they clarify; they are never the point.
4. [Topic] in BioLang — working code that applies the concept to biological data. Pipe-chained, readable, annotated.
5. Python and R Comparison — the same analysis in scipy/statsmodels and R, so you can see how the languages compare.
6. Exercises — practice problems at three difficulty levels (Foundations, Applied, Challenge).
7. Key Takeaways — the essential points to remember, in bold-and-explanation format.

Days are designed to take 1-3 hours each. Concept-heavy days (like Day 1 on distributions) are shorter. Method-heavy days (like Day 14 on logistic regression) are longer. Work at your own pace — there is no penalty for spending two days on one topic.

Prerequisites

You need:

• A computer running Windows, macOS, or Linux
• BioLang installed — see the setup section below or Appendix A
• Basic BioLang familiarity — you can write variables, use pipes, and call functions. If you have completed Practical Bioinformatics in 30 Days or the BioLang tutorials, you are ready.
• High school math — you understand addition, multiplication, fractions, and basic algebra. That is all.

You do not need:

• A statistics course (this book is the course)
• Calculus or linear algebra (we explain everything from scratch)
• Prior experience with R, Python, or any statistics software
• A powerful machine (a laptop with 4 GB of RAM handles every exercise)

If you can run bl --version and get a version number, you are ready.

The Companion Files

Every day has a companion directory with runnable code, sample data, and expected output. The structure looks like this:

biostatistics/
  days/
    day-01/
      README.md           # Day overview and instructions
      init.bl             # Setup script — run this first
      scripts/
        analysis.bl       # BioLang solutions
        analysis.py       # Python equivalent
        analysis.R        # R equivalent
      expected/
        output.txt        # Expected output for verification
      compare.md          # Side-by-side language comparison
    day-02/
      ...

To use the companion files:

1. Run init.bl first. Each day’s init script generates sample datasets, downloads reference data, or creates whatever that day’s exercises need. Run it with bl run init.bl.

2. Work through the exercises. Try to solve them yourself before looking at the solutions in scripts/.

3. Check your output. Compare your results against the files in expected/ to verify correctness. Statistical results should match within rounding tolerance.

4. Read compare.md. After completing a day in BioLang, read the comparison document to see how the same analyses look in Python and R. This is especially valuable if you plan to work in multi-language environments.

To get the companion files:

git clone https://github.com/bioras/practical-biostatistics.git
cd practical-biostatistics

Or download the ZIP from the book’s website and extract it.

Setting Up Your Environment

Full installation instructions are in Appendix A, but here is the short version:

# Install BioLang
curl -sSf https://biolang.org/install.sh | sh

# Verify it works
bl --version

# Launch the REPL to test
bl repl

On Windows, use the PowerShell installer:

irm https://biolang.org/install.ps1 | iex

If you want to run the Python comparison scripts (optional but recommended):

pip install scipy numpy pandas matplotlib statsmodels lifelines scikit-learn

If you want to run the R comparison scripts (optional but recommended):

install.packages(c("stats", "survival", "ggplot2", "dplyr", "pwr", "lme4", "boot"))

A Quick Taste

Here is what statistical analysis looks like in BioLang. This script loads gene expression data, runs a t-test between two conditions, and generates a volcano plot — all in pipe style:

# Load expression data for two conditions
let ctrl = read_csv("data/expression.csv")
let treat = read_csv("data/expression.csv")

# Run t-tests for every gene, correct for multiple testing
let results = ctrl
  |> join(treat, "gene_id")
  |> mutate("pvalue", ttest(ctrl_expr, treat_expr).p)
  |> mutate("log2fc", log2(mean(treat_expr) / mean(ctrl_expr)))
  |> mutate("padj", p_adjust(pvalue, "BH"))
  |> mutate("significant", padj < 0.05 and abs(log2fc) > 1.0)

# How many genes are differentially expressed?
results
  |> filter(significant)
  |> len()
  |> println("Differentially expressed genes: {}")

# Volcano plot
results |> volcano_plot("log2fc", "padj", "gene_id")

Twelve lines. No imports. No boilerplate. The pipe operator makes the analytical logic visible: load, join, test, correct, filter, plot. You will understand every line of this by Day 10.

Here is another example — survival analysis in three lines:

let patients = read_csv("data/clinical.csv")

patients
  |> kaplan_meier("months", "deceased", "treatment")
  |> surv_plot({title: "Overall Survival by Treatment Arm"})

And power analysis for planning your next experiment:

# Power calculation: how many samples per group?
let result = power_t_test(0.5, 0.05, 0.8)
println("Required sample size per group: {result.n}")
println("Effect size: {result.effect_size}, alpha: {result.alpha}, power: {result.power}")

BioLang’s 400+ statistical builtins mean you spend your time thinking about the biology, not fighting the syntax.

Week-by-Week Overview

Week 1: Foundations (Days 1-5)

You start where every statistical analysis starts — with the data itself. What does a distribution look like? How do you measure center and spread? What is probability, and why does it matter for hypothesis testing? Day 1 introduces distributions with real gene expression data. Day 2 covers descriptive statistics. Day 3 tackles probability and the normal distribution. Day 4 introduces sampling and the central limit theorem. Day 5 covers confidence intervals. By Friday, you have the vocabulary and intuition to understand every test that follows.

Week 2: Core Methods (Days 6-12)

Now the testing begins. Day 6 introduces hypothesis testing and p-values. Day 7 covers t-tests — one-sample, two-sample, paired — with gene expression data. Day 8 is ANOVA for comparing multiple groups. Day 9 handles non-parametric alternatives for when your data violates assumptions. Day 10 tackles chi-square and Fisher’s exact tests for categorical data. Day 11 introduces correlation and simple linear regression. Day 12 brings multiple testing correction — Bonferroni, Benjamini-Hochberg, and permutation — the single most important topic for genomics.

Week 3: Modeling (Days 13-20)

You move from testing to modeling. Day 13 covers multiple regression. Day 14 introduces logistic regression for binary outcomes. Day 15 is survival analysis — Kaplan-Meier curves and log-rank tests. Day 16 continues with Cox proportional hazards models. Day 17 introduces PCA and dimensionality reduction. Day 18 covers clustering — hierarchical, k-means, and silhouette analysis. Day 19 tackles effect sizes and confidence intervals as alternatives to p-values. Day 20 brings statistical visualization — volcano plots, Manhattan plots, Q-Q plots, forest plots, and heatmaps.

Week 4: Advanced Topics and Capstones (Days 21-27)

You tackle the hard problems. Day 21 covers bootstrap and permutation methods. Day 22 introduces Bayesian statistics with biological examples. Day 23 is power analysis and sample size calculation. Day 24 covers experimental design — randomization, blocking, batch effects. Day 25 tackles mixed models for repeated measures and nested designs. Day 26 introduces enrichment analysis — gene ontology, pathway analysis, GSEA. Day 27 covers meta-analysis for combining results across studies.

Capstone Projects (Days 28-30)

Three full projects that integrate everything you have learned. Day 28: conduct a complete RNA-seq differential expression study with quality control, normalization, testing, multiple correction, and pathway enrichment. Day 29: analyze a genome-wide association study with Manhattan plots, Q-Q plots, and genomic inflation correction. Day 30: analyze a clinical trial dataset with survival analysis, subgroup comparisons, and a statistical report suitable for publication.

Conventions Used in This Book

Throughout this book, you will see several recurring elements:

Code Blocks

BioLang code appears in fenced code blocks:

let data = [2.3, 4.1, 3.7, 5.2, 4.8]
mean(data)         # 4.02
stdev(data)        # 1.082

When a code block shows REPL interaction, lines starting with bl> are what you type:

bl> ttest([23.1, 25.4, 22.8], [19.2, 20.1, 18.7])
TTestResult { t: 4.12, df: 4, p: 0.0146 }

Shell commands use bash syntax:

bl run day07_ttest.bl

Python and R Comparisons

Multi-language comparisons appear with labeled blocks:

BioLang:

let data = read_csv("data/expression.csv")
data |> ttest(ctrl, treated) |> println()

Python:

import pandas as pd
from scipy.stats import ttest_ind
data = pd.read_csv("expression.csv")
stat, p = ttest_ind(data["ctrl"], data["treated"])
print(f"t={stat:.4f}, p={p:.4f}")

R:

data <- read.csv("expression.csv")
t.test(data$ctrl, data$treated)

Callout Boxes

Important notes, insights, and warnings appear as blockquotes throughout:

Key insight: A statistically significant result is not necessarily biologically meaningful. Always report effect sizes alongside p-values.

Clinical relevance: In oncology trials, a hazard ratio below 0.7 is typically considered clinically meaningful, regardless of the p-value.

Common pitfall: Running 20 t-tests on the same dataset without multiple testing correction gives you a 64% chance of at least one false positive at alpha = 0.05.

Exercises

Each day ends with exercises labeled by difficulty:

• Foundations — reinforce the core concept with guided problems
• Applied — use the method on a new biological dataset
• Challenge — extend the method or combine it with previous days

Key Takeaways

Each day concludes with a bulleted list of the most important points:

• The p-value is not the probability that your hypothesis is wrong. It is the probability of observing data this extreme if the null hypothesis were true. This distinction matters enormously.

A Note on the Multi-Language Approach

This book uses BioLang as its primary language because its statistical builtins let you focus on the concepts rather than the plumbing. A t-test is one function call, not a chain of imports and data manipulations. A volcano plot is one line, not thirty.

But the real world uses Python and R for most biostatistics. We include comparisons for two reasons:

1. Translation. If you already know scipy or R’s stats package, seeing the BioLang equivalent helps you learn faster. If you learn BioLang first, seeing the Python and R equivalents prepares you for collaborative work.

2. Verification. Running the same analysis in three languages and getting the same answer builds confidence. When your BioLang t-test gives p = 0.014 and your R t-test gives p = 0.014, you know you have done it right.

The compare.md file in each day’s companion directory provides a detailed side-by-side comparison. The analysis.py and analysis.R scripts are runnable equivalents you can execute and compare.

Let’s Begin

You have everything you need. The next 30 days will transform how you think about biological data — not just how to analyze it, but how to reason about uncertainty, variability, and evidence.

Day 1 starts with the most fundamental question in statistics: what does your data look like?

Turn the page. Your journey starts now.

Day 1: Why Statistics? The Story Your Data Is Trying to Tell

The Problem

In 2006, a pharmaceutical company invested over $800 million developing a promising cancer therapy. Phase II clinical trials had shown a statistically significant survival benefit in 87 patients with advanced colorectal cancer. The data looked compelling: a 38% improvement in median progression-free survival. Investors were jubilant. The company fast-tracked Phase III.

Then came the reckoning. The Phase III trial enrolled 1,200 patients across 120 medical centers. The drug performed no better than placebo. The stock price collapsed overnight. Nearly a billion dollars and eight years of research, gone — not because the science was wrong, but because the statistics were misunderstood. The Phase II “signal” was noise, amplified by a sample too small to tell the difference.

This story is not unusual. It plays out every year in laboratories, hospitals, and boardrooms around the world. Similar scenarios have unfolded with Alzheimer’s drugs, cardiovascular therapies, and anti-inflammatory agents. The difference between a breakthrough and a blunder often comes down to a few fundamental statistical concepts — concepts that any biologist, clinician, or bioinformatician can learn.

That is what this book is about. In 30 days, you will build the statistical intuition and practical skills to avoid these mistakes — whether you are designing a clinical trial, analyzing RNA-seq data, or evaluating a paper over morning coffee.

What Is Statistics?

Statistics is the science of learning from data in the presence of uncertainty. Think of it as a translator between the messy, noisy world of observations and the clean, confident conclusions you want to draw.

Imagine you are standing in a dark room, trying to understand the shape of an object by touching it with gloves on. You can feel something — ridges, curves, a rough texture — but every touch is imprecise. You might mistake a bump for an edge, or miss a hole entirely. Statistics gives you a flashlight. Not a perfect one — the beam flickers and the lens is smudged — but it is incomparably better than groping in the dark.

In biology, the “dark room” is enormous. A single human genome contains 3.2 billion base pairs. A transcriptomics experiment measures expression levels for 20,000 genes simultaneously. A clinical trial tracks hundreds of variables across thousands of patients over years. No human can intuit patterns in data this vast. Statistics provides the framework to ask precise questions and get defensible answers.

The Reproducibility Crisis

In 2012, a team at Amgen, one of the world’s largest biotechnology companies, attempted to reproduce 53 landmark studies in cancer biology — papers published in top-tier journals by respected labs. These were not obscure findings; they were studies that had shaped drug development programs and clinical practice.

They could reproduce only 6. That is an 89% failure rate.

Around the same time, Bayer HealthCare reported a similar effort. Of 67 preclinical studies they attempted to validate, roughly two-thirds could not be reproduced. The results that had guided millions of dollars in investment simply vanished when subjected to rigorous replication.

How does published, peer-reviewed science fail at this rate? The reasons are many, but they share a common root: insufficient statistical reasoning.

P-Hacking: Torturing Data Until It Confesses

One of the most insidious contributors is “p-hacking” — the practice of trying multiple analyses until one produces a statistically significant result. A researcher might:

• Test 15 different subgroups and report only the one with p < 0.05
• Remove outliers selectively until the result becomes significant
• Try multiple statistical tests and report whichever gives the smallest p-value
• Add or remove covariates until the “right” answer appears
• Decide when to stop collecting data based on whether the current result is significant

None of these practices involves fabricating data. Each individual decision might even seem reasonable in isolation. But collectively, they dramatically inflate the false positive rate. If you flip through enough combinations, you will find “significance” by pure chance — it is a mathematical certainty.

Underpowered Studies

Many published studies use sample sizes far too small to reliably detect the effects they claim to find. A study with only 12 mice per group has roughly a 20% chance of detecting a moderate treatment effect. That means 80% of real effects go undetected. But the 20% that are detected appear larger than they truly are (because only the noisiest, most extreme results cross the significance threshold), creating a distorted picture of biology.

The Garden of Forking Paths

Researchers make dozens of analytical decisions: how to clean the data, which variables to include, how to handle missing values, which test to use, whether to transform the data, how to define the outcome. Each decision is a fork in the path, and different choices lead to different results. When these choices are made after seeing the data (rather than pre-specified in an analysis plan), the researcher unconsciously navigates toward significance.

This is not an indictment of individual scientists. Most researchers receive minimal formal training in statistics. A typical biology PhD might include one semester-long course, crammed between lab rotations and qualifying exams. The result is a generation of brilliant experimentalists who treat statistical tests as black boxes — input data, output significance. Reviewers, equally uncertain about statistics, wave the paper through. The system rewards novelty over rigor.

Key insight: The reproducibility crisis is not primarily a crisis of fraud or incompetence. It is a crisis of statistical literacy. Understanding the concepts in this book is one of the most impactful things you can do for the quality of your science.

Signal vs. Noise

Here is the most fundamental question in statistics: Is the pattern I see real, or could it have happened by chance?

Consider a simple experiment. You flip a coin 10 times and get 8 heads. Is the coin biased? Your intuition says maybe — 8 out of 10 is a lot of heads. But if you do the math, a fair coin produces 8 or more heads about 5.5% of the time. That is unlikely, but not astronomically so. You might just be unlucky.

Now flip the coin 100 times and get 80 heads. Is the coin biased? Almost certainly yes. A fair coin producing 80 or more heads in 100 flips has a probability of about 0.000000000000006. You would need to flip coins continuously for the age of the universe to expect this by chance.

The pattern (80% heads) is the same in both cases. What changed is the sample size. With 10 flips, 80% heads is plausible noise. With 100 flips, 80% heads is an unmistakable signal.

This is exactly what happened with the cancer drug. In 87 patients, a 38% improvement could easily arise from random variation — which patients happened to be enrolled, how they responded to the placebo, what comorbidities they had. In 1,200 patients, the noise averages out, and the true effect (or lack thereof) becomes visible.

Common pitfall: Small studies frequently produce dramatic-looking results. This is not because small studies discover larger effects — it is because small samples are inherently noisy, and noise occasionally looks like a big signal. This phenomenon is called the “winner’s curse” and it haunts biomedical research.

The Cost of Being Wrong

In statistics, there are exactly two ways to be wrong, and they have very different consequences.

Type I Error: The False Alarm

A Type I error occurs when you conclude there is an effect when there is none. You declare the coin biased when it is actually fair. You approve a drug that does not work.

The most devastating Type I error in pharmaceutical history may be thalidomide in the 1950s. Marketed as a safe sedative for pregnant women, the drug was approved based on inadequate evidence. It caused severe birth defects in over 10,000 children worldwide. While this tragedy involved failures far beyond statistics — regulatory, ethical, and scientific — the core issue was concluding safety from data that could not support that conclusion.

A more modern example: in 2004, Merck withdrew Vioxx (rofecoxib), a blockbuster anti-inflammatory drug, after it became clear that it significantly increased heart attack risk. The drug had been on the market for five years. Post-withdrawal analysis suggested that the cardiovascular risk had been detectable in the original trial data, but was either missed or downplayed. The cost: an estimated 88,000-140,000 excess cases of heart disease in the United States alone.

Type II Error: The Missed Discovery

A Type II error occurs when you fail to detect a real effect. You declare the coin fair when it is actually biased. You reject a drug that actually works.

The canonical example is Helicobacter pylori. In 1982, Barry Marshall and Robin Warren proposed that stomach ulcers were caused by a bacterium, not by stress or spicy food as the medical establishment believed. Their initial data was compelling but their sample sizes were small. The medical community dismissed their findings for over a decade, costing millions of patients effective treatment. Marshall eventually infected himself with H. pylori, developed gastritis, and cured it with antibiotics to prove his point. He and Warren won the Nobel Prize in 2005.

Every year, real treatments are abandoned because clinical trials were too small to detect their effect. Every year, genuine biological mechanisms are dismissed because the experiment lacked statistical power. Type II errors are the silent killers of science — you never know what you missed, because the missed discovery never makes it into a journal.

How many effective cancer therapies have been shelved because the Phase II trial enrolled 40 patients instead of 400? We will never know. But the statistical tools to prevent this — power analysis and sample size calculation — are straightforward. You will learn them on Day 26.

Clinical relevance: In diagnostic testing, Type I errors produce false positives (telling a healthy person they have cancer) and Type II errors produce false negatives (telling a cancer patient they are healthy). Both are harmful, but in different ways. The balance between them is one of the central tensions in medicine.

Why Biology Needs Statistics More Than Most Fields

A physicist measuring the speed of light will get the same answer (within measurement precision) whether the experiment is run in Tokyo or Toronto, on Monday or Friday, in summer or winter. The speed of light does not have a bad day.

Biology is fundamentally different, for three reasons.

1. Biological Variability

Every organism is unique. Two genetically identical mice raised in the same cage, fed the same diet, will still differ in gene expression, tumor growth rate, immune response, and lifespan. This is not experimental error — it is the intrinsic variability of living systems. Evolution has built variability into every level of biology, from stochastic gene expression to somatic mutation to behavioral differences.

This variability means that a single measurement tells you almost nothing. If one mouse responds to a drug, you cannot conclude the drug works. If one patient’s tumor shrinks, you cannot attribute it to the treatment. You need replicates, and you need statistics to make sense of them.

2. Measurement Noise

Biological measurement is imprecise. A sequencing run introduces base-calling errors at a rate of roughly 0.1-1% per base. RNA-seq quantification depends on library preparation, read depth, alignment parameters, and normalization method. Mass spectrometry measurements fluctuate with instrument calibration, sample preparation, and ionization efficiency.

Every measurement in biology is the true signal plus some unknown amount of noise. Statistics provides the tools to separate one from the other.

3. Massive Parallel Testing

Modern biology is high-dimensional. A genome-wide association study (GWAS) tests millions of genetic variants. A differential expression analysis tests 20,000 genes. A proteomics experiment quantifies thousands of proteins. A drug screen tests hundreds of compounds.

When you test 20,000 hypotheses simultaneously, you expect 1,000 false positives by chance alone (at the conventional 5% threshold). Without proper statistical correction for multiple testing, you would drown in spurious results. This is not a theoretical concern — it is the daily reality of genomics, and getting it wrong has real consequences.

To make this concrete, consider a differential expression analysis. You measure expression of 20,000 genes in treatment versus control. Even if the treatment does absolutely nothing — affects zero genes — testing each gene at α = 0.05 will flag approximately 1,000 genes as “significant.” If you published a paper claiming these 1,000 genes are treatment-responsive, every single one would be a false positive.

The solution (multiple testing correction, which we cover on Day 15) reduces the significance threshold to account for the number of tests. In a GWAS with 1 million variants, the genome-wide significance threshold is p < 5 × 10⁻⁸ — one thousand times more stringent than the usual 0.05. Understanding why this correction is necessary, and how to apply it properly, is one of the core skills of a computational biologist.

Key insight: Biology sits at the intersection of high variability, high noise, and high dimensionality. This makes it arguably the field most in need of statistical sophistication, yet it has historically been one of the least statistically trained.

A Tour of What Lies Ahead

This book will take you from zero to practicing biostatistician in 30 days. Here is a preview of the journey:

Week 1 (Days 1-5): Foundations. You will learn to summarize data, understand distributions, reason about probability, and appreciate why sample size matters. These are the tools you need before you can test any hypothesis.

Week 2 (Days 6-10): Hypothesis Testing. You will learn confidence intervals, p-values, t-tests, and non-parametric alternatives. By Day 8, you will be able to rigorously determine whether two groups differ. By Day 10, you will know when to use (and when to avoid) parametric tests.

Week 3 (Days 11-15): Beyond Two Groups. ANOVA, chi-square tests, correlation, regression, and multiple testing correction. You will analyze multi-group experiments, categorical data, and learn why “correlation does not imply causation” is more nuanced than it sounds.

Week 4 (Days 16-20): Advanced Methods. Survival analysis, logistic regression, principal component analysis, and clustering. You will build Kaplan-Meier curves, classify patients, reduce high-dimensional data, and find natural groupings in gene expression datasets.

Week 5 (Days 21-25): Genomics Applications. Differential expression analysis, enrichment testing, multiple testing correction in practice, and Bayesian thinking. The methods that power modern computational biology.

Week 6 (Days 26-30): Real-World Practice. Power analysis and study design, meta-analysis, machine learning basics, reproducible research practices, and a capstone project that ties everything together.

By Day 8, you will know if two groups truly differ. By Day 17, you will build survival curves that predict patient outcomes. By Day 22, you will analyze differential gene expression. By Day 30, you will design a complete statistical analysis plan for a GWAS.

Each day follows the same pattern: a real-world problem that motivates the method, the conceptual framework, hands-on BioLang code, comparisons with Python and R, and exercises to cement your understanding. The emphasis throughout is on understanding — not memorizing formulas, but developing the intuition to know which method to use and why.

The Statistician’s Mindset

Before we dive into formulas and code, internalize these four questions. Ask them every time you look at data, read a paper, or plan an experiment:

1. How variable is it?

A mean without a measure of spread is almost meaningless. “Average tumor size decreased by 2 cm” sounds impressive until you learn that the standard deviation was 4 cm. Always ask: what is the spread?

2. Could chance explain this?

The human brain is wired to see patterns, even in random noise. We see faces in clouds, constellations in random stars, and trends in stock prices. Before accepting any pattern as real, quantify the probability that it arose by chance. This is the essence of hypothesis testing.

3. How big is the effect?

Statistical significance and practical significance are not the same thing. With a large enough sample, you can detect a difference of 0.001 grams in tumor weight with p < 0.001. But is a one-milligram difference clinically meaningful? Always report effect sizes alongside p-values.

4. Is my sample representative?

If you study the genetics of heart disease using only patients from a single hospital in Boston, your results may not generalize to patients in rural India. If you select only the “best” cell lines for your experiment, your conclusions may not extend to primary cells. Sampling bias is the silent assassin of biomedical research.

Putting the Mindset into Practice

These four questions are not abstract philosophy. They are a practical checklist:

You will encounter these questions again and again throughout this book. By Day 30, they will be second nature — the automatic mental checklist of a statistically literate scientist.

Key insight: Statistics is not a set of tests to run after the experiment. It is a way of thinking that should inform every stage — from study design to data collection to analysis to interpretation. The best time to consult a statistician is before you collect a single data point.

The Burden of Proof

In everyday life, we make decisions based on intuition, anecdote, and authority. “My grandmother smoked until 95, so smoking cannot be that bad.” “This supplement worked for my friend, so it must be effective.” “The famous professor says this treatment works, so it must.”

Science demands a higher standard. The burden of proof rests on the claimant. If you claim a drug works, you must demonstrate it with evidence strong enough to withstand scrutiny. If you claim a gene is associated with a disease, you must show that the association is unlikely to be a coincidence.

Statistics provides the machinery for this burden of proof. It forces you to be explicit about your assumptions, quantify your uncertainty, and acknowledge the limits of your data. It is, in essence, formalized humility.

Consider the claim “Vitamin D supplements reduce cancer risk.” An anecdote is worthless: your uncle took vitamin D and did not get cancer. A small observational study is weak: 50 people who took vitamin D had fewer cancers than 50 who did not — but maybe the vitamin D group was healthier to begin with (confounding). A large randomized controlled trial with 25,000 participants, pre-registered outcomes, and proper statistical analysis is strong evidence. Each step up the ladder requires more statistical sophistication.

The hierarchy of evidence is, fundamentally, a hierarchy of statistical rigor:

This book will equip you to evaluate and produce evidence at every level of this hierarchy.

The Numbers Tell a Story

To bring this all together, let us look at a real-world scenario that illustrates every concept from today.

A research group publishes a paper claiming that a new biomarker predicts response to immunotherapy. Their study: 24 patients (12 responders, 12 non-responders). They measure the biomarker level in each patient and find a “statistically significant” difference (p = 0.03).

Here is what a statistical thinker would ask:

How variable is it? The biomarker levels range from 2 to 200 ng/mL. The standard deviation within each group is enormous — nearly as large as the difference between groups. The signal is weak relative to the noise.

Could chance explain it? With only 12 per group and high variability, the p-value of 0.03 is fragile. If you removed two extreme patients, it becomes 0.12. The result is not robust.

How big is the effect? The difference in medians is 15 ng/mL, but the overlap between groups is substantial. Many responders have lower biomarker levels than many non-responders. The effect size (Cohen’s d) is only 0.4 — a “small to medium” effect.

Is the sample representative? All patients came from a single institution, were predominantly male, and had a specific tumor subtype. Whether the biomarker works in a broader population is unknown.

A naive reader sees “p < 0.05, significant.” A statistically literate reader sees a fragile, underpowered result from a non-representative sample with a modest effect size. These are different conclusions from the same data.

A Preview of the Tools

Throughout this book, you will use BioLang to perform statistical analyses. Here is a tiny glimpse of what Day 2 will look like — just to whet your appetite:

# Tomorrow, you'll summarize 10,000 quality scores in one line:
# let stats = summary(quality_scores)
#
# And visualize them instantly:
# histogram(quality_scores, {bins: 50, title: "Sequencing Quality Distribution"})

But today is about the why, not the how. The tools are only as good as the thinking behind them. A researcher who understands why a t-test exists will use it correctly even with imperfect software. A researcher who merely knows how to call a t-test function will misuse it regularly, regardless of how elegant the software is.

Exercises

Exercise 1: The Newspaper Test

Find a news article reporting a scientific or medical finding (e.g., “Coffee reduces cancer risk by 15%”). Write down your answers to these four questions:

• (a) What was the sample size? If the article does not mention it, what does that tell you?
• (b) Could the result be due to chance? What would you need to know to answer this?
• (c) Is the effect size meaningful in practice? A 2% reduction in cancer risk sounds different from a 50% reduction.
• (d) Is the sample representative of the population you care about? Who was studied, and who was not?

If the article does not provide enough information to answer these questions, that itself is informative. Most science journalism omits sample sizes, effect sizes, and confidence intervals — precisely the information you need to evaluate the claim.

Exercise 2: Coin Flip Thought Experiment

Without doing any math, estimate the following:

• If you flip a fair coin 20 times, what is the probability of getting exactly 10 heads?
• What about 15 or more heads out of 20?
• What about 20 heads in a row?

Write down your guesses. We will revisit this on Day 4 with the tools to compute exact answers, and you can see how well your intuition calibrated.

Exercise 3: Reproducibility Reflection

Think about a result from your own work (or a paper you have read) that you found surprising or striking. List three reasons why the result might fail to reproduce if someone repeated the experiment. For each reason, identify whether it is:

• Statistical (sample size, random variation, multiple testing)
• Methodological (different protocols, reagent lots, equipment)
• Biological (different cell lines, patient populations, environmental conditions)

Exercise 4: Type I vs Type II in Your Field

Identify one example each of a Type I error (false positive) and a Type II error (false negative) that would be particularly damaging in your area of biology. For each:

• Describe the scenario concretely
• Estimate the consequences (financial, clinical, scientific)
• State which error type you consider more dangerous in your context, and why

Exercise 5: Spotting P-Hacking

A paper reports testing a drug on patients across 8 different cancer subtypes. Only one subtype shows a significant result (p = 0.04). The paper’s title highlights this positive finding. What statistical concerns should this raise? How many false positives would you expect by chance when testing 8 subtypes at α = 0.05?

Key Takeaways

• Statistics is the science of learning from data in the presence of uncertainty — it is essential, not optional, for biological research.
• The reproducibility crisis is largely a statistical literacy crisis: most landmark findings fail replication due to underpowered studies, p-hacking, and misunderstood tests.
• Signal vs. noise is the fundamental statistical question: the same percentage difference can be meaningful or meaningless depending on sample size.
• Type I errors (false positives) waste resources and can cause harm; Type II errors (false negatives) cause missed discoveries. Neither can be eliminated — only managed.
• Biology is uniquely challenging for statistics due to inherent biological variability, measurement noise, and massive parallel testing.
• The statistician’s mindset asks four questions: How variable? Could chance explain it? How big is the effect? Is the sample representative?
• Statistics should inform every stage of research, from design through interpretation — not just the analysis phase.

What’s Next

Tomorrow, we roll up our sleeves and meet real data. You will learn to summarize 10,000 numbers into a handful of meaningful statistics — means, medians, standard deviations, and more. You will discover why the mean is a liar when outliers are present, why box plots reveal truths that histograms hide, and how a single command in BioLang can tell you whether a sequencing run is worth analyzing or should be thrown away. Day 2 is where the hands-on work begins.

Day 2: Your Data at a Glance — Descriptive Statistics

The Problem

Dr. Sarah Chen stares at her screen. The Illumina NovaSeq 6000 finished its run overnight, and now she has 10,247 quality scores — one for each tile on the flow cell. Her PI needs a decision by the morning meeting: is this data usable, or do they need to re-run the library, burning another $4,000 and two days?

She cannot read 10,247 numbers. She cannot scroll through them and develop an intuition. She has five minutes before the meeting starts. What she needs is a way to compress 10,247 numbers into a handful of meaningful summaries that answer three questions: What is the typical quality? How much does it vary? Are there any red flags?

This is the job of descriptive statistics. They are the first thing you compute, every time, before you run a single test. Get them wrong — or skip them — and everything downstream is built on sand.

What Are Descriptive Statistics?

Descriptive statistics are summaries. Think of them as a movie trailer for your data. The full movie (the raw dataset) might be two hours long, but the trailer gives you the genre, the tone, and the key plot points in two minutes. A good trailer does not lie about the movie. A good set of descriptive statistics does not lie about the data.

There are three things you need to know about any dataset:

1. Center — Where is the “middle” of the data?
2. Spread — How far do values range from that center?
3. Shape — Is the data symmetric? Skewed? Are there outliers?

Why Do These Three Matter?

Center tells you what is “normal.” If the mean quality score of your sequencing run is Q35, you know tiles are performing well. If it is Q15, something went wrong. Without center, you have no reference point — you cannot say whether a single observation is typical or alarming. In clinical genomics, center defines the baseline: what is the typical variant allele frequency in this cohort? What is the average coverage across your target regions?

Spread tells you how much you can trust the center. Two experiments might both report a mean IC50 of 12 nM, but one has values ranging from 11 to 13 (tight, reproducible) while the other ranges from 2 to 45 (noisy, unreliable). The center is the same — the spread tells you completely different stories. High spread means your next measurement could land anywhere; low spread means you can make confident predictions. In RNA-seq, high within-group variance buries real differential expression in noise.

Shape tells you which statistical tools are safe to use. Almost every standard statistical test (t-test, ANOVA, linear regression) assumes the data is approximately bell-shaped (normally distributed). If your data is heavily right-skewed — which gene expression nearly always is — those tests give wrong answers. Shape also reveals hidden structure: a bimodal distribution might mean you have two distinct populations mixed together (e.g., responders and non-responders to a drug). Ignoring shape is the single most common reason published biomedical results fail to replicate.

Key insight: Center without spread is meaningless (“the average temperature is 72°F” tells you nothing if the range is -20 to 160). Spread without shape is dangerous (a standard deviation assumes symmetry, but skewed data violates this). You always need all three.

Every statistical analysis starts here. If you skip descriptive statistics and jump straight to hypothesis testing, you are performing surgery without an examination.

Measures of Center

Mean (Arithmetic Average)

The mean is the balance point. If your data values were weights placed along a ruler, the mean is the spot where the ruler would balance perfectly.

Formula: x̄ = (1/n) ∑ xᵢ

The mean uses every data point, which is both its strength and its weakness. It is the most efficient estimator of center when data is symmetric with no outliers. But it is exquisitely sensitive to extreme values.

Example: Five gene expression values (FPKM): 12, 15, 14, 13, 16. Mean = 14.0. Reasonable.

Now add one highly expressed gene: 12, 15, 14, 13, 16, 5000. Mean = 845.0. The mean has been dragged from 14 to 845 by a single outlier. It no longer represents “typical” expression.

Common pitfall: In genomics, gene expression distributions are heavily right-skewed. Reporting mean FPKM/TPM values without acknowledging this skew is misleading. The median is almost always a better summary for expression data.

Median (Middle Value)

The median is the value that splits the data in half: 50% of observations fall below it, 50% above. Sort the data and pick the middle number (or average the two middle numbers if n is even).

For our outlier-contaminated expression data: sorted = 12, 13, 14, 15, 16, 5000. Median = (14 + 15) / 2 = 14.5. The outlier barely matters.

The median is robust — it resists the pull of extreme values. This makes it the preferred measure of center for skewed distributions, which are the norm in biology.

Mode (Most Frequent Value)

The mode is the most common value. It is most useful for categorical or discrete data: the most common blood type, the most frequent variant allele, the peak of a histogram.

For continuous data, the mode is the peak of the density curve. Bimodal distributions (two peaks) arise in biology more often than you might expect — for instance, CpG methylation levels often cluster near 0% and 100%, reflecting fully unmethylated and fully methylated states.

When to Use Each

Key insight: Always report both mean and median. If they differ substantially, your data is skewed, and the median is the more honest summary.

Measures of Spread

Knowing the center is not enough. Two datasets can have identical means and wildly different behaviors. Consider drug response in two patient cohorts — both might have a mean survival of 12 months, but in one cohort everyone lives 11-13 months, while in the other, half die in 2 months and half live 22 months. The clinical implications are completely different.

Range

The simplest measure: maximum minus minimum. It tells you the total extent of the data but is completely determined by the two most extreme points.

Variance and Standard Deviation

Variance measures the average squared distance from the mean:

Variance: s² = (1 / (n-1)) ∑ (xᵢ - x̄)²

Standard Deviation: s = √(s²)

We divide by (n-1) rather than n (Bessel’s correction) because a sample underestimates the true population variance. The standard deviation is in the same units as the data, making it more interpretable than variance.

Rule of thumb: For normally distributed data, about 68% of values fall within 1 SD of the mean, 95% within 2 SDs, and 99.7% within 3 SDs. If a quality score is more than 3 SDs below the mean, something is wrong with that tile.

Interquartile Range (IQR)

The IQR is the range of the middle 50% of the data: Q3 - Q1, where Q1 is the 25th percentile and Q3 is the 75th percentile.

Like the median, the IQR is robust to outliers. It is the foundation of the box plot and a standard measure of spread for skewed data.

Coefficient of Variation (CV)

CV = (SD / Mean) × 100%. The CV expresses variability relative to the mean, making it useful for comparing spread across measurements with different scales.

Example: If gene A has mean expression 1000 TPM with SD 100, and gene B has mean 10 TPM with SD 5, which is more variable? Gene A has higher SD, but gene B has a higher CV (50% vs 10%). Gene B’s expression is relatively more variable.

Shape: Skewness and Kurtosis

Skewness

Skewness measures asymmetry. A skewness of 0 means perfectly symmetric (like the normal distribution). Positive skewness means a right tail (most values cluster low, with a few very high values). Negative skewness means a left tail.

Biological reality: Most biological measurements are positively skewed. Gene expression, protein abundance, cell counts, read depths — all tend to have many low values and a few very high ones. This is because multiplicative processes (gene regulation cascades, exponential growth) naturally produce right-skewed distributions.

Kurtosis

Kurtosis measures the “tailedness” of the distribution — how likely extreme values are compared to a normal distribution. High kurtosis means heavy tails (more outliers than expected). Low kurtosis means light tails.

In genomics, variant allele frequency distributions often have high kurtosis, reflecting the mixture of common variants (clustered near 50%) and rare variants (near 0%).

Key insight: If skewness is far from 0 or kurtosis is far from 3, your data is not normally distributed, and parametric tests that assume normality may give misleading results. We will explore this deeply on Day 3.

Descriptive Statistics in BioLang

Let us return to Dr. Chen’s problem. She has 10,247 quality scores and five minutes.

Loading and Summarizing Data

set_seed(42)
# Simulate sequencing quality scores (Phred scale, typically 0-40)
let quality_scores = rnorm(10247, 32.5, 3.8)

# One-line comprehensive summary
let stats = summary(quality_scores)
print(stats)
# Output:
#   n:        10247
#   mean:     32.48
#   median:   32.51
#   sd:       3.81
#   min:      17.23
#   max:      47.12
#   q1:       29.93
#   q3:       35.07
#   iqr:      5.14
#   skewness: -0.01
#   kurtosis: 2.99

That single call answers Dr. Chen’s three questions immediately. Mean quality is 32.5 (good — above the Phred 30 threshold). The SD is 3.8 (moderate spread). Skewness is near zero (symmetric). No red flags. The run is usable.

Computing Individual Statistics

# Individual measures of center
let avg = mean(quality_scores)       # 32.48
let med = median(quality_scores)     # 32.51
let mod = mode(quality_scores)       # 32 (rounded to nearest integer)

# Individual measures of spread
let sd = stdev(quality_scores)       # 3.81
let v = variance(quality_scores)     # 14.52
let r = range_stat(quality_scores)   # [17.23, 47.12]
let q = quantile(quality_scores, [0.25, 0.5, 0.75])  # [29.93, 32.51, 35.07]
let i = iqr(quality_scores)          # 5.14

# Shape
let sk = skewness(quality_scores)    # -0.01
let ku = kurtosis(quality_scores)    # 2.99

print("Mean: {avg}, Median: {med}")
print("SD: {sd}, IQR: {i}")
print("Skewness: {sk}, Kurtosis: {ku}")

The summary() Function

For a more detailed report:

let report = summary(quality_scores)
print(report)
# Output:
#   Variable:    quality_scores
#   Count:       10247
#   Mean:        32.48
#   Std Dev:     3.81
#   Min:         17.23
#   25%:         29.93
#   50%:         32.51
#   75%:         35.07
#   Max:         47.12
#   Skewness:    -0.01
#   Kurtosis:    2.99
#   CV:          11.73%
#   SE(mean):    0.038

Working with Real Sequencing Data

set_seed(42)
# In practice, you'd load from a file:
# let scores = read_csv("data/expression.csv") |> column("phred_score")

# Simulate a problematic run with bimodal quality
let good_tiles = rnorm(8000, 33.0, 2.5)
let bad_tiles = rnorm(2247, 18.0, 4.0)
let mixed_scores = good_tiles + bad_tiles

let mixed_stats = summary(mixed_scores)
print(mixed_stats)
# Mean: 29.7 — looks okay at first glance
# Median: 32.1 — the median reveals the truth: most tiles are fine
# Skewness: -1.4 — strongly left-skewed, warning sign!

# The mean alone would have hidden the problem.
# Always look at the full distribution.

Common pitfall: Relying on the mean alone can hide bimodal distributions. The “bad tile” problem above is common in sequencing — a mean of 29.7 looks passable, but 22% of tiles are failing. Always visualize.

Visualization: Always Plot Before You Test

Numbers tell part of the story. Plots tell the rest. Anscombe’s Quartet — four datasets with identical means, variances, and correlations but wildly different structures — demonstrates why you must always look at your data.

Histograms

A histogram bins your data and counts frequencies. It reveals the distribution’s shape at a glance.

# Basic histogram (uses quality_scores from block above — click "Run All Above + This")
histogram(quality_scores, {bins: 50, title: "Sequencing Quality Distribution"})

# Histogram of the problematic run (uses mixed_scores from earlier block)
histogram(mixed_scores, {bins: 50, title: "Bimodal Quality — Problem Run"})
# This immediately reveals the two peaks that the mean masked.

The good run produces a single symmetric peak around 32-33. The problem run shows two distinct peaks — one centered at 18 and one at 33. No summary statistic captures this as effectively as the histogram.

Box Plots

A box plot displays the median (center line), IQR (box), and whiskers (typically 1.5 × IQR). Points beyond the whiskers are marked as individual outliers.

set_seed(42)
# Compare quality across lanes
let lane1 = rnorm(2000, 33.0, 2.0)
let lane2 = rnorm(2000, 31.5, 3.5)
let lane3 = rnorm(2000, 28.0, 5.0)

let bp_table = table({"Lane 1": lane1, "Lane 2": lane2, "Lane 3": lane3})
boxplot(bp_table, {title: "Quality Score Distribution by Lane"})
# Lane 3 immediately stands out: lower median, wider spread, more outliers.

Box plots shine when comparing groups. You can see at a glance which lane is problematic, without computing a single number.

Combining Plots and Stats

set_seed(42)
# Full QC report in a few lines
let scores = rnorm(10000, 32.5, 3.8)

# Side-by-side: histogram + box plot
histogram(scores, {bins: 40, title: "Quality Score Distribution"})
boxplot(scores, {title: "Quality Score Box Plot"})

# Summary table
let stats = summary(scores)
print("QC Decision: {if stats.mean > 30 then 'PASS' else 'FAIL'}")
print("Mean Q-score: {stats.mean:.1}")
print("Tiles below Q20: {scores |> filter(|s| s < 20) |> len()}")

Python and R Equivalents

For those coming from Python or R, here are the equivalent operations.

Python:

import numpy as np
from scipy import stats

scores = np.random.normal(32.5, 3.8, 10247)

# Measures of center
np.mean(scores)          # 32.48
np.median(scores)        # 32.51
stats.mode(scores)       # most frequent

# Measures of spread
np.std(scores, ddof=1)   # 3.81 (ddof=1 for sample SD)
np.var(scores, ddof=1)   # 14.52
np.percentile(scores, [25, 50, 75])
stats.iqr(scores)        # 5.14

# Shape
stats.skew(scores)       # -0.01
stats.kurtosis(scores)   # -0.01 (scipy uses excess kurtosis)

# Comprehensive summary
import pandas as pd
pd.Series(scores).describe()

R:

scores <- rnorm(10247, mean = 32.5, sd = 3.8)

# Measures of center
mean(scores)       # 32.48
median(scores)     # 32.51

# Measures of spread
sd(scores)         # 3.81
var(scores)        # 14.52
quantile(scores, c(0.25, 0.5, 0.75))
IQR(scores)        # 5.14

# Shape (requires moments package)
library(moments)
skewness(scores)   # -0.01
kurtosis(scores)   # 2.99

# Comprehensive summary
summary(scores)

# Visualization
hist(scores, breaks = 50, main = "Quality Distribution")
boxplot(scores, main = "Quality Box Plot")

Worked Example: The QC Decision

Let us walk through Dr. Chen’s complete analysis.

set_seed(42)
# Step 1: Load the data
let scores = rnorm(10247, 32.5, 3.8)

# Step 2: Quick summary
let stats = summary(scores)

# Step 3: QC criteria
let pass_threshold = 30.0     # Minimum acceptable mean quality
let fail_fraction_limit = 0.10 # Max 10% tiles below Q20

# Step 4: Compute QC metrics
let tiles_below_q20 = scores |> filter(|s| s < 20) |> len()
let fraction_below_q20 = tiles_below_q20 / len(scores)

# Step 5: Decision
let qc_pass = stats.mean >= pass_threshold and fraction_below_q20 <= fail_fraction_limit

print("=== Sequencing Run QC Report ===")
print("Total tiles:          {len(scores)}")
print("Mean quality:         {stats.mean:.2}")
print("Median quality:       {stats.median:.2}")
print("Std deviation:        {stats.sd:.2}")
print("Tiles below Q20:      {tiles_below_q20} ({fraction_below_q20 * 100:.1}%)")
print("QC Result:            {if qc_pass then 'PASS' else 'FAIL'}")
print("================================")

# Step 6: Visualize
histogram(scores, {bins: 50, title: "Run QC: Quality Score Distribution"})

This takes about 10 seconds to run. Dr. Chen has her answer well before the meeting.

Exercises

Exercise 1: Gene Expression Summary

Compute descriptive statistics for a set of gene expression values and identify the best measure of center.

# Gene expression values (TPM) for 20 genes
let expression = [0.5, 1.2, 3.4, 5.1, 8.7, 12.3, 15.0, 22.4, 45.6, 78.9,
                  120.5, 250.3, 0.1, 2.8, 6.5, 0.3, 1100.0, 33.2, 0.8, 18.5]

# TODO: Compute mean, median, stdev, skewness
# TODO: Which measure of center best represents "typical" expression? Why?
# TODO: Create a histogram. What shape do you see?

Exercise 2: Comparing Sequencing Runs

Two sequencing runs produced quality scores. Determine which run is better.

set_seed(42)
let run_a = rnorm(5000, 31.0, 2.5)
let run_b = rnorm(5000, 33.0, 6.0)

# TODO: Compute summary() for both runs
# TODO: Which has higher mean? Which has lower variability?
# TODO: Compute the CV for each. Which is more consistent?
# TODO: Create side-by-side box plots
# TODO: If you could only pick one run, which would you choose and why?

Exercise 3: Outlier Detection

Identify outliers using the IQR method (values below Q1 - 1.5IQR or above Q3 + 1.5IQR).

# Protein abundance measurements with some suspicious values
let protein = [45.2, 48.1, 50.3, 47.8, 49.5, 46.7, 51.2, 48.9,
               150.0, 47.3, 49.8, 46.1, 50.5, 48.4, 3.0, 49.1]

# TODO: Compute Q1, Q3, IQR
# TODO: Calculate lower and upper fences
# TODO: Identify which values are outliers
# TODO: Compute mean with and without outliers — how much does it change?

Exercise 4: Multi-Sample QC Dashboard

Build a QC summary for multiple samples.

set_seed(42)
let samples = {
    "Sample_A": rnorm(1000, 34.0, 2.0),
    "Sample_B": rnorm(1000, 29.0, 4.0),
    "Sample_C": rnorm(1000, 32.0, 3.0),
    "Sample_D": rnorm(1000, 20.0, 5.0)
}

# TODO: Loop through samples, compute summary() for each
# TODO: Flag any sample with mean < 25 or CV > 20%
# TODO: Create a box plot comparing all four samples

Key Takeaways

• Descriptive statistics compress large datasets into interpretable summaries of center, spread, and shape.
• The mean is efficient but outlier-sensitive; the median is robust. For skewed biological data, prefer the median.
• Standard deviation measures absolute spread; IQR is robust to outliers; CV enables comparison across different scales.
• Skewness and kurtosis reveal whether your data’s shape matches the assumptions of common statistical tests.
• Always visualize your data with histograms and box plots before computing any test. Summary statistics can hide multimodal distributions, outliers, and other structural features.
• summary() in BioLang provides comprehensive descriptive statistics in a single function call.
• Descriptive statistics are not optional preliminaries — they are the foundation of every analysis.

What’s Next

You now know how to summarize data, but you may have noticed something: we keep saying “normally distributed” and “skewed” without precisely defining what a distribution is. Tomorrow, on Day 3, we dive into the mathematical shapes that biological data follows — the normal distribution, the log-normal, the Poisson, and the binomial. You will learn why gene expression data refuses to be normal, why read counts follow a Poisson process (sort of), and how to test whether your data fits the distribution you think it does. Understanding distributions is the key to choosing the right statistical test — and avoiding the wrong one.

Day 3: Distributions — The Shape of Biological Variation

The Problem

Dr. James Park has RNA-seq data from 12 tumor samples. He needs to identify genes that are differentially expressed between treatment and control groups. He knows the standard approach: run a t-test on each gene. The t-test assumes data is normally distributed, so he checks his data.

Gene FPKM values range from 0 to 50,000. Most genes sit near zero, with a handful expressed at astronomical levels. The histogram looks nothing like a bell curve — it is a massive spike at zero with a long right tail stretching to infinity. He runs the t-test anyway. Out of 20,000 genes, 4,700 come back “significant” at p < 0.05. That is 23.5% of all genes, far more than seems biologically plausible for this modest treatment.

Something is deeply wrong. The t-test assumed his data was symmetric and bell-shaped. It was neither. The test produced garbage because the assumption was violated. If Dr. Park had understood distributions — the subject of today’s chapter — he would have known to transform his data or use a different test entirely.

What Is a Distribution?

A distribution is a recipe that tells you how likely each possible value is. Think of it as a map of a city, where the height at each point represents how many people live there. Downtown is a tall spike (many people). The suburbs are a gentle slope. The surrounding farmland is nearly flat.

Every dataset has an underlying distribution — the theoretical shape that generated the data you observe. When you draw a histogram of your data, you are estimating this shape from a finite sample. With 10 data points, the histogram is choppy and unreliable. With 10,000, it smooths out and begins to reveal the true underlying curve.

Why does this matter? Because every statistical test makes assumptions about the distribution of your data. The t-test assumes normality. The chi-square test assumes expected counts are large enough. The Poisson regression assumes counts follow a Poisson process. Violate the assumption, and the test’s guarantees evaporate.

Key insight: A statistical test is a contract. It says: “If your data follows distribution X, then I guarantee my conclusions are reliable with probability Y.” Break the contract, and you get no guarantees.

The Normal Distribution

The normal distribution — the bell curve — is the most famous distribution in statistics, and for good reason. It arises naturally whenever many small, independent effects add together.

Properties

The normal distribution is defined by two parameters:

• μ (mu): the mean, which determines the center
• σ (sigma): the standard deviation, which determines the width

The curve is perfectly symmetric around μ. It extends infinitely in both directions, though values far from the mean are exceedingly unlikely.

The 68-95-99.7 Rule

If a measurement falls more than 3 standard deviations from the mean, it is either a genuine outlier or something went wrong.

Biological Examples of Normality

The normal distribution is a good model for:

• Measurement error: Technical replicates of the same sample tend to be normally distributed.
• Height and weight in a homogeneous population (though mixtures of populations are not normal).
• Blood pressure readings in healthy adults.
• Quantitative traits influenced by many genes of small effect (additive genetic model).
• Log-transformed gene expression (more on this below).

When Data Is NOT Normal

The normal distribution is a terrible model for:

• Raw gene expression (FPKM, TPM, counts) — heavily right-skewed
• Read counts — discrete, non-negative, often zero-inflated
• Allele frequencies — bounded between 0 and 1
• Survival times — always positive, typically right-skewed
• Any data with a hard boundary (concentrations cannot be negative)

set_seed(42)
# Generate and visualize normal data
let heights = rnorm(5000, 170, 8)

histogram(heights, {bins: 50, title: "Adult Heights (cm) — Normal Distribution"})
let stats = summary(heights)
print("Mean: {stats.mean:.1}, Median: {stats.median:.1}, Skewness: {stats.skewness:.3}")
# Mean and median nearly identical; skewness near zero — hallmarks of normality

The Log-Normal Distribution

If gene expression is not normal, what is it? In most cases, it is log-normal: the data itself is skewed, but the logarithm of the data is normally distributed.

Why Gene Expression Is Log-Normal

Gene regulation is a cascade of multiplicative processes. A transcription factor binds (or doesn’t), an enhancer activates (fold-change), mRNA is stabilized (half-life multiplied), ribosomes translate at varying rates (multiplied efficiency). When effects multiply rather than add, the result is log-normal, not normal.

This is a mathematical fact: if X = Y₁ × Y₂ × … × Yₙ and the Y values are independent, then log(X) = log(Y₁) + log(Y₂) + … + log(Yₙ). Sums of independent variables tend toward normal (by the Central Limit Theorem), so log(X) is approximately normal, meaning X is log-normal.

The Log-Transform Trick

This is why bioinformaticians routinely log-transform expression data before analysis:

set_seed(42)
# Simulate gene expression (log-normal)
let log_expr = rnorm(5000, 3.0, 2.0)
let expression = log_expr |> map(|x| 2.0 ** x)  # 2^x to simulate FPKM

# Raw expression: heavily skewed
histogram(expression, {bins: 50, title: "Raw FPKM — Right Skewed"})
let raw_stats = summary(expression)
print("Raw — Mean: {raw_stats.mean:.1}, Median: {raw_stats.median:.1}, Skew: {raw_stats.skewness:.2}")

# Log2-transformed: approximately normal
let log2_expr = expression |> map(|x| log2(x + 1))  # +1 to handle zeros
histogram(log2_expr, {bins: 50, title: "log2(FPKM+1) — Approximately Normal"})
let log_stats = summary(log2_expr)
print("Log2 — Mean: {log_stats.mean:.1}, Median: {log_stats.median:.1}, Skew: {log_stats.skewness:.2}")

After log-transformation, the mean and median converge, skewness drops toward zero, and the histogram looks bell-shaped. Now parametric tests are appropriate.

Clinical relevance: Differential expression tools like DESeq2 and edgeR work with counts and model them with the negative binomial distribution, but many downstream analyses (clustering, PCA, visualization) require log-transformed data. Understanding why is essential for correct analysis.

The Poisson Distribution

The Poisson distribution models the number of events that occur in a fixed interval, when events happen independently at a constant average rate.

Properties

• Parameter: λ (lambda) — the average rate
• Support: 0, 1, 2, 3, … (non-negative integers)
• Mean = Variance = λ (this is the key property)

Biological Examples

The Overdispersion Problem

In theory, RNA-seq counts should be Poisson. In practice, biological replicates show more variability than Poisson predicts — the variance exceeds the mean. This is overdispersion, caused by biological variability between samples.

The solution is the negative binomial distribution, which adds a dispersion parameter to allow variance > mean. This is why DESeq2 uses negative binomial, not Poisson.

set_seed(42)
# Poisson distribution for mutation counts

# Average 3.5 mutations per megabase in a tumor
let mutation_counts = rpois(1000, 3.5)

histogram(mutation_counts, {bins: 15, title: "Mutations per Megabase (Poisson, lambda=3.5)"})

# Verify mean ~ variance (Poisson property)
print("Mean: {mean(mutation_counts):.2}")
print("Variance: {variance(mutation_counts):.2}")
# Both should be close to 3.5

# Probability of seeing 10+ mutations in a region (hypermutation?)
let p_hyper = 1.0 - ppois(9, 3.5)
print("P(10+ mutations): {p_hyper:.4}")
# Very low — a region with 10+ mutations is genuinely unusual

The Binomial Distribution

The binomial distribution models the number of successes in a fixed number of independent trials, each with the same probability of success.

Properties

• Parameters: n (number of trials), p (probability of success)
• Support: 0, 1, 2, …, n
• Mean = np, Variance = np(1-p)

Biological Examples

Hardy-Weinberg Equilibrium

For a biallelic locus with allele frequencies p and q = 1-p, Hardy-Weinberg predicts genotype frequencies:

• AA: p²
• Aa: 2pq
• aa: q²

Deviations from HWE can indicate selection, population structure, or genotyping error.

set_seed(42)
# Genotype frequencies under Hardy-Weinberg
let p = 0.3  # frequency of allele A
let q = 1.0 - p

let freq_AA = p * p         # 0.09
let freq_Aa = 2.0 * p * q   # 0.42
let freq_aa = q * q         # 0.49

print("Expected genotype frequencies:")
print("  AA: {freq_AA:.3}")
print("  Aa: {freq_Aa:.3}")
print("  aa: {freq_aa:.3}")

# Simulate genotyping 500 individuals
let n_individuals = 500
let AA_count = rbinom(1, n_individuals, freq_AA) |> sum()

# Probability of observing exactly k heterozygotes
let k = 200
let p_exact = dbinom(k, n_individuals, freq_Aa)
print("P(exactly {k} heterozygotes in {n_individuals}): {p_exact:.6}")

# Probability of 230+ heterozygotes (possible HWE violation?)
let p_excess = 1.0 - pbinom(229, n_individuals, freq_Aa)
print("P(230+ heterozygotes): {p_excess:.4}")

Checking Your Distribution: Diagnostic Tools

Before running any parametric test, verify that your data matches its assumptions. Here are the essential tools.

Q-Q Plots

A Q-Q (quantile-quantile) plot compares your data’s quantiles against the theoretical quantiles of a reference distribution (usually normal). If data follows the reference distribution, points fall on a straight diagonal line. Deviations reveal departures from the assumed shape.

set_seed(42)
# Q-Q plot for normal data (should be a straight line)
let normal_data = rnorm(500, 0, 1)
qq_plot(normal_data, {title: "Q-Q Plot: Normal Data"})

# Q-Q plot for log-normal data (curved — not normal!)
let lognormal_data = normal_data |> map(|x| exp(x))
qq_plot(lognormal_data, {title: "Q-Q Plot: Log-Normal Data (Raw)"})

# Q-Q plot after log-transform (straight again)
let transformed = lognormal_data |> map(|x| log(x))
qq_plot(transformed, {title: "Q-Q Plot: Log-Normal Data (After Log Transform)"})

Reading a Q-Q plot:

• Points on the line: data matches the assumed distribution
• Upward curve at the right end: right skew (heavy right tail)
• Downward curve at the left end: left skew (heavy left tail)
• S-shape: heavy tails on both sides (high kurtosis)

Shapiro-Wilk Test

The Shapiro-Wilk test formally tests whether data is normally distributed. A small p-value (< 0.05) means the data is significantly non-normal.

set_seed(42)
# Test normality visually with Q-Q plots and histograms
let normal_data = rnorm(200, 50, 10)
let skewed_data = rnorm(200, 3, 1) |> map(|x| exp(x))

# For normal data: Q-Q plot should show points on the diagonal
qq_plot(normal_data, {title: "Q-Q: Normal Data"})
let stats_normal = summary(normal_data)
print("Normal data — Skewness: {stats_normal.skewness:.4}")
# Skewness near 0: consistent with normality

# For skewed data: Q-Q plot will curve away from the diagonal
qq_plot(skewed_data, {title: "Q-Q: Skewed Data"})
let stats_skewed = summary(skewed_data)
print("Skewed data — Skewness: {stats_skewed.skewness:.4}")
# High skewness: definitely not normal

Common pitfall: The Shapiro-Wilk test is very sensitive with large samples. With n = 10,000, it will reject normality for data that is “close enough” to normal for practical purposes. For large samples, rely more on Q-Q plots and skewness/kurtosis values than on the Shapiro-Wilk p-value.

Histogram with Density Overlay

Overlay the theoretical density curve on your histogram to visually assess fit:

set_seed(42)
let data = rnorm(2000, 100, 15)

# Histogram with normal density overlay
histogram(data, {bins: 40, title: "Data with Normal Fit Overlay"})
density(data, {title: "Kernel Density Estimate"})

Distribution Summary Table

Python and R Equivalents

Python:

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# Normal distribution
x = np.linspace(-4, 4, 1000)
plt.plot(x, stats.norm.pdf(x, 0, 1))

# Poisson
counts = np.random.poisson(lam=3.5, size=1000)

# Binomial
genotypes = np.random.binomial(n=500, p=0.42, size=1000)

# Q-Q plot
stats.probplot(data, dist="norm", plot=plt)

# Shapiro-Wilk test
stat, p = stats.shapiro(data)

# Distribution fitting
params = stats.norm.fit(data)  # MLE fit

R:

# Normal distribution
x <- seq(-4, 4, length.out = 1000)
plot(x, dnorm(x, 0, 1), type = "l")

# Poisson
counts <- rpois(1000, lambda = 3.5)

# Binomial
genotypes <- rbinom(1000, size = 500, prob = 0.42)

# Q-Q plot
qqnorm(data)
qqline(data)

# Shapiro-Wilk test
shapiro.test(data)

# Density overlay
hist(data, freq = FALSE)
curve(dnorm(x, mean(data), sd(data)), add = TRUE, col = "red")

Why This Matters for Testing

Here is the critical connection between today’s material and the rest of the book:

Choosing the wrong test because you assumed the wrong distribution is one of the most common errors in computational biology. Dr. Park’s mistake from our opening scenario — running a t-test on raw FPKM values — is committed daily in bioinformatics labs around the world.

Key insight: The distribution is not a detail. It is the foundation. Get it right, and your downstream analysis is trustworthy. Get it wrong, and no amount of sophisticated testing can rescue your conclusions.

Exercises

Exercise 1: Identify the Distribution

For each dataset, determine the most appropriate distribution and justify your choice.

set_seed(42)
# Dataset A: Sequencing read counts per gene
let dataset_a = rpois(1000, 25)

# Dataset B: Patient blood pressure
let dataset_b = rnorm(500, 120, 15)

# Dataset C: Gene expression (raw TPM)
let log_vals = rnorm(2000, 2.0, 1.5)
let dataset_c = log_vals |> map(|x| 10.0 ** x)

# TODO: For each dataset, create a histogram and Q-Q plot
# TODO: Check normality visually with qq_plot() and histogram()
# TODO: For dataset C, try log-transforming and re-check
# TODO: State which distribution best describes each and why

Exercise 2: The Poisson Check

Verify whether mutation count data follows a Poisson distribution by checking the mean-variance relationship.

set_seed(42)
# Scenario 1: Pure Poisson (technical replicates)
let technical = rpois(500, 8.0)

# Scenario 2: Overdispersed (biological replicates)
# Simulate by mixing Poisson with varying lambda
let lambdas = rnorm(500, 8.0, 3.0) |> map(|x| max(x, 0.1))

# TODO: Compute mean and variance for technical replicates
# TODO: Are they approximately equal? (Poisson property)
# TODO: For overdispersed data, compute the dispersion ratio (variance/mean)
# TODO: What does a ratio >> 1 tell you about the data?

Exercise 3: Transform and Test

Take a skewed dataset, find the right transformation, and verify normality.

set_seed(42)
let protein_abundance = rnorm(300, 4, 1.2) |> map(|x| exp(x))

# TODO: Plot histogram of raw data
# TODO: Check normality with qq_plot() — is it normal?
# TODO: Apply log transform
# TODO: Plot histogram of transformed data
# TODO: Check normality of transformed data with qq_plot()
# TODO: Compare skewness before and after

Exercise 4: Distribution Detective

A collaborator gives you mystery data. Identify its distribution using all tools from today.

set_seed(42)
# Mystery data — what distribution is this?
let mystery = rbinom(1000, 50, 0.15)

# TODO: Compute summary()
# TODO: Create histogram
# TODO: Try Q-Q plot against normal
# TODO: Note: the data is discrete. What distributions produce discrete data?
# TODO: Estimate the parameters and identify the distribution

Key Takeaways

• A distribution is the theoretical shape describing how likely each value is. Every dataset has one, and every statistical test assumes one.
• The normal distribution arises from additive effects and is defined by mean and standard deviation. It is appropriate for measurement error and many physiological traits.
• Gene expression is NOT normal — it is log-normal because gene regulation is multiplicative. Always log-transform before using parametric tests.
• The Poisson distribution models count data (reads, mutations) with the key property that mean equals variance. When variance exceeds the mean (overdispersion), use the negative binomial instead.
• The binomial distribution models fixed trials with a success probability — relevant for genotype frequencies and allele sampling.
• Q-Q plots are the most informative visual diagnostic for distribution checking. The Shapiro-Wilk test provides a formal hypothesis test for normality.
• Choosing the right distribution is not optional — it determines which statistical tests are valid and which will produce misleading results.

What’s Next

You now understand the shapes that biological data takes. But distributions describe what values are likely — which is just another way of saying they describe probabilities. Tomorrow, on Day 4, we formalize probability itself. You will learn to compute the chance that a child inherits a BRCA1 mutation, understand why a positive genetic test might mean less than you think (Bayes’ theorem will surprise you), and discover why the prosecutor’s fallacy has sent innocent people to prison. Probability is the language of uncertainty, and uncertainty is the native tongue of biology.

Day 4: Probability — Quantifying Uncertainty

The Problem

Maria and David sit in a genetic counselor’s office. The air is still. Maria has just learned she carries a pathogenic BRCA1 mutation — a variant that dramatically increases lifetime risk of breast and ovarian cancer. They are planning to start a family and they need answers.

What is the probability their child inherits the mutation? If the child inherits it, what is the probability she develops breast cancer by age 70? They are considering preimplantation genetic testing — if the test says the embryo is mutation-free, how confident can they be? The counselor pulls out a notepad and begins writing probabilities.

This is not an abstract exercise. These numbers will determine whether Maria and David proceed with natural conception, pursue IVF with genetic screening, or consider adoption. The difference between a 50% risk and a 5% risk changes lives. Understanding how to compute, combine, and interpret probabilities is not just statistics — in genetics, it is clinical care.

What Is Probability?

Probability is a number between 0 and 1 that quantifies how likely an event is to occur. Zero means impossible. One means certain. Everything interesting happens in between.

Think of probability as a weather forecast. When the forecaster says “70% chance of rain,” she means: in historical situations with similar atmospheric conditions, it rained about 70% of the time. She is not saying it will rain exactly 70% of the total rainfall. She is quantifying uncertainty about a future event using past data and models.

In biology, we use probability constantly:

• The probability a child inherits a specific allele from a heterozygous parent: 0.5
• The probability a random human carries at least one pathogenic BRCA1 variant: roughly 1/400
• The probability that a drug produces a response in a given cancer type: varies, but measured in clinical trials
• The probability that a sequencing read is mapped incorrectly: the mapping quality score encodes this directly

Basic Rules of Probability

The Addition Rule

The probability of event A or event B occurring depends on whether they can both happen at the same time.

Mutually exclusive events (cannot co-occur): P(A or B) = P(A) + P(B)

A person’s blood type is A, B, AB, or O. These are mutually exclusive — you cannot be both type A and type B simultaneously. So P(A or B) = P(A) + P(B) = 0.40 + 0.11 = 0.51.

Non-mutually exclusive events (can co-occur): P(A or B) = P(A) + P(B) - P(A and B)

A patient might have diabetes, hypertension, or both. If P(diabetes) = 0.10, P(hypertension) = 0.30, and P(both) = 0.05, then P(either) = 0.10 + 0.30 - 0.05 = 0.35. You subtract the overlap to avoid double-counting.

The Multiplication Rule

The probability of event A and event B both occurring depends on whether they are independent.

Independent events (one does not affect the other): P(A and B) = P(A) × P(B)

The probability that two unrelated people both carry a BRCA1 mutation: P(carrier) × P(carrier) = (1/400) × (1/400) = 1/160,000.

Dependent events (one affects the other): P(A and B) = P(A) × P(B|A)

where P(B|A) is the probability of B given that A has occurred. If a mother carries BRCA1 (event A), the probability her daughter both inherits it (B) and develops cancer by 70 (C) is: P(B and C) = P(B|A) × P(C|B) = 0.5 × 0.72 = 0.36.

The Complement Rule

P(not A) = 1 - P(A)

The probability of not inheriting the mutation is 1 - 0.5 = 0.5. The probability a diagnostic test does not give a false positive is 1 - (false positive rate). Simple but powerful — often easier to compute the complement and subtract.

Key insight: Most probability errors in biology come from confusing independent and dependent events, or from forgetting to subtract the overlap in non-mutually exclusive events. Write out the formula before plugging in numbers.

Conditional Probability

Conditional probability is the probability of an event given that another event has occurred. It is written P(A|B) and read “the probability of A given B.”

This is where most people’s intuition breaks down, because P(A|B) is not the same as P(B|A).

The Critical Distinction

• P(cancer | BRCA1 mutation) ≈ 0.72 — If you carry BRCA1, your lifetime breast cancer risk is about 72%.
• P(BRCA1 mutation | cancer) ≈ 0.05 — If you have breast cancer, the probability it is due to BRCA1 is only about 5%.

These are completely different numbers answering completely different questions. Confusing them is called the inverse probability fallacy, and it has real consequences in medicine, law, and genetics.

The Prosecutor’s Fallacy

In forensic genetics, the prosecutor’s fallacy works like this: “The probability of this DNA match occurring by chance is 1 in 10 million. Therefore, the probability the defendant is innocent is 1 in 10 million.”

This is logically wrong. P(match | innocent) is not P(innocent | match). In a city of 8 million people, you would expect roughly one other person to match by chance. If the only evidence is the DNA match, the probability of innocence might be closer to 50%, not 1 in 10 million.

The same fallacy appears in genetic testing: “The test is 99% accurate” does not mean a positive result is 99% likely to be correct. The answer depends on how common the condition is — which brings us to Bayes’ theorem.

Common pitfall: Never interpret P(result | hypothesis) as P(hypothesis | result). This mistake is ubiquitous in clinical genetics, drug development, and forensics. Bayes’ theorem is the correction.

Bayes’ Theorem

Bayes’ theorem provides the mathematical machinery to flip conditional probabilities. Given P(B|A), it computes P(A|B):

P(A|B) = P(B|A) × P(A) / P(B)

In words: the probability of A given B equals the probability of B given A, times the prior probability of A, divided by the total probability of B.

The Diagnostic Test Example

This is the single most important application of Bayes’ theorem in biomedical science. Work through it carefully.

Setup:

• A genetic disease has a prevalence of 1 in 1,000 (P(disease) = 0.001)
• A diagnostic test has 99% sensitivity: P(positive | disease) = 0.99
• The test has 95% specificity: P(negative | no disease) = 0.95, so P(positive | no disease) = 0.05

Question: A patient tests positive. What is the probability they actually have the disease?

Intuition says: The test is 99% sensitive and 95% specific — surely a positive result means ~95-99% chance of disease?

Bayes says:

P(disease | positive) = P(positive | disease) × P(disease) / P(positive)

P(positive) = P(positive | disease) × P(disease) + P(positive | no disease) × P(no disease) P(positive) = 0.99 × 0.001 + 0.05 × 0.999 P(positive) = 0.00099 + 0.04995 = 0.05094

P(disease | positive) = 0.00099 / 0.05094 = 0.0194

A positive test result means only a 1.94% chance of actually having the disease.

How can a “99% accurate” test give such a low positive predictive value? Because the disease is rare. In 100,000 people, 100 have the disease (99 test positive) and 99,900 are healthy (4,995 test positive by mistake). Of the 5,094 total positives, only 99 are true positives.

PPV = 99 / 5,094 = 1.94%. The false positives overwhelm the true positives.

Clinical relevance: This is why population-wide genetic screening for rare conditions produces mostly false positives. It is also why a confirmatory test with a different methodology is always required. Understanding PPV is essential for anyone interpreting genetic test results.

Bayes in BioLang

# Diagnostic test calculator using Bayes' theorem
let prevalence = 0.001      # P(disease) = 1 in 1,000
let sensitivity = 0.99      # P(positive | disease)
let specificity = 0.95      # P(negative | no disease)
let fpr = 1.0 - specificity # P(positive | no disease) = 0.05

# Total probability of testing positive
let p_positive = sensitivity * prevalence + fpr * (1.0 - prevalence)

# Positive Predictive Value (PPV)
let ppv = (sensitivity * prevalence) / p_positive

# Negative Predictive Value (NPV)
let p_negative = (1.0 - sensitivity) * prevalence + specificity * (1.0 - prevalence)
let npv = (specificity * (1.0 - prevalence)) / p_negative

print("=== Diagnostic Test Analysis ===")
print("Prevalence:    {prevalence}")
print("Sensitivity:   {sensitivity}")
print("Specificity:   {specificity}")
print("P(positive):   {p_positive:.4}")
print("PPV:           {ppv:.4} ({ppv * 100:.1}%)")
print("NPV:           {npv:.6} ({npv * 100:.4}%)")
print("")
print("Interpretation: A positive result means only a {ppv * 100:.1}% chance of disease.")
print("A negative result means a {npv * 100:.4}% chance of being disease-free.")

How Prevalence Changes Everything

# Show how PPV changes with disease prevalence
let sensitivity = 0.99
let specificity = 0.95

let prevalences = [0.0001, 0.001, 0.01, 0.05, 0.10, 0.20, 0.50]

print("Prevalence | PPV")
print("-----------|--------")
for prev in prevalences {
    let fpr = 1.0 - specificity
    let p_pos = sensitivity * prev + fpr * (1.0 - prev)
    let ppv = (sensitivity * prev) / p_pos
    print("  {prev:.4}   | {ppv * 100:.1}%")
}
# At 0.01% prevalence: PPV = 0.2% (nearly all positives are false)
# At 50% prevalence: PPV = 95.2% (most positives are true)

This table is one of the most important results in clinical statistics. It explains why screening tests work well for common conditions but poorly for rare ones.

Probability Distributions Revisited

On Day 3, we met distributions as shapes. Now we can understand them as probability functions.

Discrete Distributions

For discrete random variables (counts, genotypes), the probability mass function (PMF) gives P(X = k) — the probability of observing exactly the value k.

# Binomial PMF: probability of exactly k successes in n trials
# Scenario: 4 children, mother is BRCA1 carrier (p = 0.5)
let n_children = 4
let p_inherit = 0.5

print("Number of children inheriting BRCA1:")
for k in 0..5 {
    let prob = dbinom(k, n_children, p_inherit)
    print("  {k} children: {prob:.4} ({prob * 100:.1}%)")
}
# 0: 6.25%, 1: 25.0%, 2: 37.5%, 3: 25.0%, 4: 6.25%

Cumulative Distribution

The cumulative distribution function (CDF) gives P(X ≤ k) — the probability of observing k or fewer.

# What's the probability that at most 1 of 4 children inherits the mutation?
let p_at_most_1 = pbinom(1, 4, 0.5)
print("P(0 or 1 child inherits): {p_at_most_1:.4}")  # 0.3125

# What's the probability at least 1 inherits?
let p_at_least_1 = 1.0 - pbinom(0, 4, 0.5)
print("P(at least 1 inherits): {p_at_least_1:.4}")  # 0.9375

Continuous Distributions

For continuous random variables (gene expression, blood pressure), the probability density function (PDF) does not give P(X = k) (which is always zero for continuous variables). Instead, probabilities are computed over intervals using the CDF.

# Blood pressure: normal with mean 120, SD 15
let mu = 120.0
let sigma = 15.0

# P(BP > 140) — hypertension threshold
let p_hypertension = 1.0 - pnorm(140, mu, sigma)
print("P(BP > 140): {p_hypertension:.4}")  # ~0.0912

# P(100 < BP < 130)
let p_normal_range = pnorm(130, mu, sigma) - pnorm(100, mu, sigma)
print("P(100 < BP < 130): {p_normal_range:.4}")

# What BP value has only 5% of people above it?
let bp_95th = qnorm(0.95, mu, sigma)
print("95th percentile BP: {bp_95th:.1}")  # ~144.7

Hardy-Weinberg as a Probability Model

Hardy-Weinberg equilibrium (HWE) is one of the most elegant applications of probability in genetics. For a biallelic locus with allele frequencies p (allele A) and q = 1-p (allele a), random mating produces genotypes with probabilities:

This model treats each allele transmission as an independent random event — like drawing two alleles from a bag with replacement.

# Test for Hardy-Weinberg equilibrium
# Observed genotype counts from a population study
let observed_AA = 45
let observed_Aa = 210
let observed_aa = 245
let total = observed_AA + observed_Aa + observed_aa  # 500

# Estimate allele frequencies
let p_A = (2.0 * observed_AA + observed_Aa) / (2.0 * total)
let p_a = 1.0 - p_A
print("Estimated allele frequencies: p(A) = {p_A:.3}, p(a) = {p_a:.3}")

# Expected counts under HWE
let expected_AA = p_A * p_A * total
let expected_Aa = 2.0 * p_A * p_a * total
let expected_aa = p_a * p_a * total
print("Expected: AA={expected_AA:.1}, Aa={expected_Aa:.1}, aa={expected_aa:.1}")
print("Observed: AA={observed_AA}, Aa={observed_Aa}, aa={observed_aa}")

# Chi-square test for HWE (we'll cover this formally on Day 13)
let chi2 = (observed_AA - expected_AA) ** 2 / expected_AA +
           (observed_Aa - expected_Aa) ** 2 / expected_Aa +
           (observed_aa - expected_aa) ** 2 / expected_aa
print("Chi-square statistic: {chi2:.3}")
print("Deviation from HWE: {if chi2 > 3.84 then 'Significant' else 'Not significant'}")

Carrier Probability Calculations

Returning to Maria and David’s consultation:

# Genetic counseling probability calculator

# Maria is a BRCA1 carrier (heterozygous)
# David is not a carrier (assumed)

# Autosomal dominant: 50% chance each child inherits
let p_inherit = 0.5

# They want 3 children
let n_children = 3

print("=== Genetic Counseling: BRCA1 Inheritance ===")
print("")

# Probability none of 3 children inherit
let p_none = dbinom(0, n_children, p_inherit)
print("P(no children inherit): {p_none:.4} ({p_none * 100:.1}%)")

# Probability exactly 1 inherits
let p_one = dbinom(1, n_children, p_inherit)
print("P(exactly 1 inherits):  {p_one:.4} ({p_one * 100:.1}%)")

# Probability at least 1 inherits
let p_at_least_one = 1.0 - p_none
print("P(at least 1 inherits): {p_at_least_one:.4} ({p_at_least_one * 100:.1}%)")

print("")
print("=== Conditional Cancer Risk ===")
# If a daughter inherits BRCA1, lifetime breast cancer risk ~ 72%
let p_cancer_given_brca = 0.72

# P(inherits AND develops cancer)
let p_inherit_and_cancer = p_inherit * p_cancer_given_brca
print("P(daughter inherits AND gets cancer): {p_inherit_and_cancer:.3} ({p_inherit_and_cancer * 100:.1}%)")

# P(daughter gets cancer by 70 | she is female)
# Must account for 50% chance of being female
let p_affected_daughter = 0.5 * p_inherit * p_cancer_given_brca
print("P(random child is affected daughter): {p_affected_daughter:.3} ({p_affected_daughter * 100:.1}%)")

Python and R Equivalents

Python:

from scipy import stats
import numpy as np

# Binomial probabilities
stats.binom.pmf(k=2, n=4, p=0.5)     # P(X = 2)
stats.binom.cdf(k=1, n=4, p=0.5)     # P(X <= 1)

# Normal probabilities
stats.norm.cdf(140, loc=120, scale=15) # P(X <= 140)
stats.norm.ppf(0.95, loc=120, scale=15) # 95th percentile

# Poisson
stats.poisson.pmf(k=5, mu=3.5)        # P(X = 5)
stats.poisson.cdf(k=9, mu=3.5)        # P(X <= 9)

# Bayes calculation
prevalence = 0.001
sensitivity = 0.99
fpr = 0.05
p_pos = sensitivity * prevalence + fpr * (1 - prevalence)
ppv = (sensitivity * prevalence) / p_pos

R:

# Binomial
dbinom(2, size = 4, prob = 0.5)    # P(X = 2)
pbinom(1, size = 4, prob = 0.5)    # P(X <= 1)

# Normal
pnorm(140, mean = 120, sd = 15)    # P(X <= 140)
qnorm(0.95, mean = 120, sd = 15)   # 95th percentile

# Poisson
dpois(5, lambda = 3.5)             # P(X = 5)
ppois(9, lambda = 3.5)             # P(X <= 9)

# Bayes
prevalence <- 0.001
sensitivity <- 0.99
fpr <- 0.05
p_pos <- sensitivity * prevalence + fpr * (1 - prevalence)
ppv <- (sensitivity * prevalence) / p_pos

Exercises

Exercise 1: The Prenatal Test

A prenatal screening test for Down syndrome has sensitivity 95% and specificity 97%. The prevalence of Down syndrome is approximately 1 in 700 live births for a 30-year-old mother.

# TODO: Compute the PPV — if the test is positive, what is the true probability?
# TODO: Compute the NPV — if the test is negative, how reassuring is it?
# TODO: How does the PPV change for a 40-year-old mother (prevalence ~1 in 100)?
let prevalence_30 = 1.0 / 700.0
let prevalence_40 = 1.0 / 100.0
let sensitivity = 0.95
let specificity = 0.97

Exercise 2: Multiple Independent Events

A patient takes three independent diagnostic tests for a condition. Each test has 90% sensitivity.

# TODO: What is P(all three tests are positive | disease)?
# TODO: What is P(at least one test is negative | disease)?
# TODO: What is P(all three tests are negative | disease)?
# TODO: If the patient tests positive on all three, and prevalence is 1%,
#       what is the updated probability they have the disease?
let sensitivity = 0.90
let specificity = 0.95
let prevalence = 0.01

Exercise 3: Carrier Frequency

Cystic fibrosis is autosomal recessive with a carrier frequency of approximately 1 in 25 among Europeans.

let carrier_freq = 1.0 / 25.0

# TODO: What is P(both parents are carriers)?
# TODO: If both are carriers, P(affected child)?
# TODO: P(random European couple has an affected child)?
# TODO: If one parent is a confirmed carrier and the other is untested,
#       what is P(their child is affected)?

Exercise 4: Sequencing Error

A variant caller reports a SNV at a position covered by 50 reads. 5 reads support the variant allele.

# Is this a real variant or sequencing error?
# Assume sequencing error rate = 0.01 per base per read

let n_reads = 50
let n_alt = 5
let error_rate = 0.01

# TODO: Under the null (all errors), what is P(5+ alt reads)?
# Use the binomial: P(X >= 5 | n=50, p=0.01) = 1 - pbinom(4, 50, 0.01)
# TODO: Is this consistent with error alone, or likely a real variant?

Key Takeaways

• Probability quantifies uncertainty on a 0-to-1 scale. It is the mathematical language of statistics.
• The addition rule handles “or” questions; the multiplication rule handles “and” questions. Whether events are mutually exclusive or independent changes the formula.
• Conditional probability P(A|B) is NOT the same as P(B|A). Confusing them is the source of the prosecutor’s fallacy and misinterpretation of diagnostic tests.
• Bayes’ theorem is the bridge from P(B|A) to P(A|B). It shows that a positive test for a rare disease usually means the patient is healthy — the positive predictive value depends critically on prevalence.
• In genetics, probability models like Hardy-Weinberg equilibrium and Mendelian inheritance translate directly into binomial probability calculations.
• Always compute both PPV and NPV when interpreting diagnostic or screening tests. Sensitivity and specificity alone are insufficient.

What’s Next

Today you learned to compute probabilities for individual events. But in real experiments, you do not observe single events — you observe samples drawn from populations. Tomorrow, on Day 5, we tackle the crucial question of sampling: Why does sample size matter so much? You will see the Central Limit Theorem in action — watching the distribution of sample means magically approach normality even when the underlying data is wildly skewed. You will understand what “statistical power” means and why a study with 20 patients per arm is almost guaranteed to fail. Day 5 is the bridge between probability theory and the practical reality of experimental design.

Day 5: Sampling, Bias, and Why n Matters

The Problem

Dr. Elena Vasquez is the lead biostatistician for a pharmaceutical company. A new immunotherapy drug has shown promising results in cell lines and mouse models. Now the clinical team is designing the Phase II trial and they want her sign-off on the sample size.

The clinical lead proposes 20 patients per arm — treatment and placebo. “It’s faster, cheaper, and we can get to Phase III sooner,” he argues. Elena runs the numbers and shakes her head. With 20 patients per arm and the expected effect size, the trial has only a 23% chance of detecting the drug’s benefit even if it truly works. That means a 77% chance of concluding the drug is ineffective when it actually saves lives.

She recommends 200 patients per arm. The clinical lead winces at the cost — $12 million more and 18 extra months of enrollment. But Elena is firm: “Would you rather spend $12 million now and know the answer, or spend $50 million on a Phase III that was doomed from the start because Phase II was too small to see the signal?”

This tension — between the cost of collecting more data and the cost of drawing wrong conclusions from too little — is the central drama of experimental design. Today you will understand why sample size is not a bureaucratic detail but the most consequential decision in any study.

What Are Populations and Samples?

The Population

The population is the complete set of items you want to understand. It is usually too large, too expensive, or physically impossible to measure in its entirety.

The Sample

The sample is the subset you actually observe. Everything you learn comes from the sample, but everything you want to know is about the population. Statistics is the bridge between the two.

The quality of the bridge depends entirely on two factors:

1. How the sample was selected (bias)
2. How large the sample is (precision)

Key insight: A large biased sample is worse than a small unbiased one. The 1936 Literary Digest poll surveyed 2.4 million people and predicted Alf Landon would win the presidential election in a landslide. George Gallup surveyed 50,000 and correctly predicted Roosevelt. The Literary Digest sample was drawn from telephone directories and automobile registrations — overrepresenting wealthy voters. Size could not compensate for bias.

The Sampling Distribution

This is one of the most important concepts in all of statistics, and it is the one that most students find counterintuitive at first.

Imagine you draw a sample of 30 patients, measure their blood pressure, and compute the mean. You get 125 mmHg. Now imagine you draw a different sample of 30 patients and compute the mean. You might get 121 mmHg. A third sample: 128 mmHg.

If you repeated this process thousands of times — each time drawing 30 patients and computing the mean — you would get thousands of sample means. These sample means form a distribution called the sampling distribution of the mean.

The sampling distribution is NOT the distribution of individual data points. It is the distribution of a statistic (like the mean) computed from repeated samples.

Seeing It in Action

set_seed(42)
# Simulate the sampling distribution

# The "population": blood pressure values (slightly right-skewed)
let population = rnorm(100000, 125, 18)

# Draw 1000 samples of size 30, compute mean of each
let sample_means_30 = []
for i in 0..1000 {
    let s = sample(population, 30)
    sample_means_30 = sample_means_30 + [mean(s)]
}

# Draw 1000 samples of size 200, compute mean of each
let sample_means_200 = []
for i in 0..1000 {
    let s = sample(population, 200)
    sample_means_200 = sample_means_200 + [mean(s)]
}

# Compare the distributions
print("Population:         mean = {mean(population):.1}, SD = {stdev(population):.1}")
print("Sample means (n=30):  mean = {mean(sample_means_30):.1}, SD = {stdev(sample_means_30):.1}")
print("Sample means (n=200): mean = {mean(sample_means_200):.1}, SD = {stdev(sample_means_200):.1}")

histogram(sample_means_30, {bins: 40, title: "Sampling Distribution (n=30)"})
histogram(sample_means_200, {bins: 40, title: "Sampling Distribution (n=200)"})

Two crucial observations:

1. Both sampling distributions are centered at the true population mean (~125). Samples are unbiased estimators.
2. The n=200 distribution is much narrower than n=30. Larger samples give more precise estimates.

The Central Limit Theorem

The Central Limit Theorem (CLT) is perhaps the single most important result in statistics. It says:

Regardless of the shape of the population distribution, the sampling distribution of the mean approaches a normal distribution as sample size increases.

This is remarkable. The underlying data can be skewed, bimodal, uniform, or any shape at all. As long as you take large enough samples and compute means, those means will be approximately normally distributed.

Demonstrating the CLT

set_seed(42)
# Create a wildly non-normal population: exponential (very right-skewed)
let skewed_pop = rnorm(100000, 2, 1) |> map(|x| exp(x))

# The population is extremely skewed
histogram(skewed_pop, {bins: 50, title: "Population: Exponential (Very Skewed)"})
let pop_stats = summary(skewed_pop)
print("Population skewness: {pop_stats.skewness:.2}")

# Sample means with n=5 (still somewhat skewed)
let means_n5 = []
for i in 0..2000 {
    let s = sample(skewed_pop, 5)
    means_n5 = means_n5 + [mean(s)]
}
histogram(means_n5, {bins: 50, title: "Sample Means, n=5"})
print("n=5 skewness: {skewness(means_n5):.2}")

# Sample means with n=30 (approaching normal)
let means_n30 = []
for i in 0..2000 {
    let s = sample(skewed_pop, 30)
    means_n30 = means_n30 + [mean(s)]
}
histogram(means_n30, {bins: 50, title: "Sample Means, n=30"})
print("n=30 skewness: {skewness(means_n30):.2}")

# Sample means with n=100 (very close to normal)
let means_n100 = []
for i in 0..2000 {
    let s = sample(skewed_pop, 100)
    means_n100 = means_n100 + [mean(s)]
}
histogram(means_n100, {bins: 50, title: "Sample Means, n=100"})
print("n=100 skewness: {skewness(means_n100):.2}")

# Verify normality visually with Q-Q plot
qq_plot(means_n100, {title: "Q-Q Plot: Sample Means n=100"})

Watch the skewness drop toward zero as n increases. By n=100, the sampling distribution is indistinguishable from a normal curve, even though the underlying data is wildly skewed.

Key insight: The CLT is why the normal distribution dominates statistics. Even when individual observations are non-normal, means of samples are approximately normal. Since most statistical tests are fundamentally about comparing means, the normal distribution is the right reference distribution for the test statistic — even when the raw data is not normal.

When Does the CLT “Kick In”?

The speed of convergence to normality depends on how non-normal the population is:

The “n ≥ 30” rule of thumb is a rough guideline, not a universal truth.

Standard Error: The Precision of Your Estimate

The standard deviation of the sampling distribution has a special name: the standard error (SE).

SE = SD / √n

This formula encodes the fundamental relationship between sample size and precision:

• Double your sample size -> SE decreases by a factor of sqrt(2) ~ 1.41
• Quadruple your sample size -> SE halves
• To cut SE in half, you need 4 times as many observations

set_seed(42)
# Demonstrate how SE shrinks with sample size
let population_sd = 18.0  # Blood pressure SD

let sample_sizes = [5, 10, 20, 30, 50, 100, 200, 500, 1000]

print("Sample Size | Theoretical SE | Observed SE")
print("------------|----------------|------------")

for n in sample_sizes {
    let theoretical_se = population_sd / sqrt(n)

    # Simulate to verify
    let means = []
    for i in 0..1000 {
        let s = rnorm(n, 125, population_sd)
        means = means + [mean(s)]
    }
    let observed_se = stdev(means)

    print("  {n:>6}    |     {theoretical_se:>6.2}     |    {observed_se:.2}")
}

With 20 patients, your estimate of mean blood pressure could easily be off by 8 mmHg — enough to misclassify a treatment as effective or ineffective. With 200 patients, you are unlikely to be off by more than 2.5 mmHg.

Common pitfall: Researchers often confuse SD and SE. The SD describes variability among individual observations. The SE describes precision of the sample mean. They answer different questions. Report the right one. SD for describing data; SE for describing the precision of an estimate.

Types of Bias

Sample size controls precision, but even infinite precision cannot fix a biased sample. Bias is a systematic error that pushes your estimate in a consistent direction.

Selection Bias

Your sample is not representative of the population you want to study.

Example: A study of gene expression in breast cancer recruits patients only from a single academic medical center. These patients tend to have more advanced disease (referral bias), are more likely to be white (geographic bias), and have better follow-up (compliance bias). The results may not generalize to community hospitals or diverse populations.

Genomics example: If you study “healthy controls” by recruiting university employees, your sample overrepresents educated, relatively affluent people — not the general population.

Survivorship Bias

You only observe subjects who “survived” some selection process, missing those who did not.

Classic example: During WWII, the military examined bullet holes in returning planes and planned to add armor where holes were most common. Statistician Abraham Wald pointed out the error: they were only seeing planes that survived. The areas with no holes were where planes had been hit and crashed. Armor should go where holes were absent.

Biological example: If you study long-term cancer survivors to find prognostic biomarkers, you miss the patients who died quickly. Your biomarkers will predict survival among survivors, not among all patients.

Ascertainment Bias

The way you identify subjects systematically skews who gets included.

Example: A study finds that children with autism have more genetic variants than controls. But the autistic children were ascertained through clinical evaluation (which involves deep phenotyping and genetic testing), while controls were population-based. The ascertainment process itself led to more thorough variant discovery in cases.

Measurement Bias

The way you measure introduces systematic error.

Example: A technician consistently reads gel bands as slightly brighter than they are. All expression measurements are systematically inflated. If this bias is constant across all samples, relative comparisons are still valid. If it varies between conditions (e.g., the technician knows which samples are treatment), it corrupts everything.

Genomics example: GC content bias in sequencing — regions with extreme GC content are systematically under-represented in coverage, biasing any analysis that depends on read depth.

Publication Bias

Studies with significant results are more likely to be published than studies with null results. The published literature systematically overestimates effect sizes.

Example: 20 groups test whether gene X is associated with disease Y. One group (by chance) finds p < 0.05 and publishes. The other 19 find nothing and file the results away. The published literature now says gene X is associated with disease Y, but the full evidence says otherwise.

Key insight: Bias cannot be fixed by increasing sample size. A biased study with 10,000 subjects gives you a very precise wrong answer. Always evaluate bias before interpreting results.

Why n Matters: The Power Preview

Statistical power is the probability of detecting a real effect when it exists. It is 1 minus the Type II error rate (1 - β). Convention targets 80% power, meaning a 20% chance of missing a real effect.

Power depends on four factors:

1. Effect size — How large is the true difference? Bigger effects are easier to detect.
2. Sample size (n) — More data = more power.
3. Variability (σ) — Less noise = more power.
4. Significance threshold (α) — More stringent threshold = less power.

Think of detecting a treatment effect as hearing a whisper in a crowd. The whisper is the signal (effect size). The crowd noise is variability. Adding more listeners (larger n) helps. Making the crowd quieter (reducing variability) helps. Demanding absolute certainty before you will believe you heard something (lower α) makes it harder.

Simulating Power

set_seed(42)
# Simulate a clinical trial to understand power

# True effect: treatment mean = 125 (control = 130, lower is better)
let control_mean = 130.0
let treatment_mean = 125.0  # 5 mmHg real difference
let sd = 18.0

# Function to run one trial and check if we detect the difference
# Returns 1 if p < 0.05, 0 otherwise
let run_trial = fn(n_per_arm) {
    let control = rnorm(n_per_arm, control_mean, sd)
    let treatment = rnorm(n_per_arm, treatment_mean, sd)
    let result = ttest(treatment, control)
    if result.p_value < 0.05 { 1 } else { 0 }
}

# Run 1000 simulated trials for different sample sizes
let sizes = [20, 50, 100, 200, 500]
print("n per arm | Estimated Power")
print("----------|----------------")

for n in sizes {
    let detections = 0
    for i in 0..1000 {
        detections = detections + run_trial(n)
    }
    let power = detections / 1000.0
    print("  {n:>5}   |     {power * 100:.1}%")
}

Typical results:

This is Dr. Vasquez’s argument in numbers. With 20 patients per arm, the trial has a 77% chance of producing a false negative — concluding the drug does not work when it does. That is not an experiment; it is a waste of money.

The Bootstrap: Estimation Without Formulas

The bootstrap is a resampling method that estimates the sampling distribution empirically. Instead of relying on mathematical formulas, you resample your data with replacement thousands of times and compute your statistic each time.

The bootstrap is invaluable when:

• The formula for the standard error is unknown or complicated
• The CLT may not apply (small n, skewed data)
• You want confidence intervals for any statistic (median, correlation, ratio)

set_seed(42)
# Bootstrap estimation of the standard error of the median

# Original sample: 50 gene expression values
let expression = rnorm(50, 3.0, 1.5) |> map(|x| exp(x))

let observed_median = median(expression)
print("Observed median: {observed_median:.2}")

# Bootstrap: resample with replacement 5000 times
let boot_medians = []
for i in 0..5000 {
    let resample = sample(expression, len(expression))
    boot_medians = boot_medians + [median(resample)]
}

# Bootstrap SE
let boot_se = stdev(boot_medians)
print("Bootstrap SE of median: {boot_se:.2}")

# Bootstrap 95% confidence interval (percentile method)
let ci_lower = quantile(boot_medians, 0.025)
let ci_upper = quantile(boot_medians, 0.975)
print("95% Bootstrap CI: [{ci_lower:.2}, {ci_upper:.2}]")

histogram(boot_medians, {bins: 50, title: "Bootstrap Distribution of the Median"})

Key insight: The bootstrap treats your sample as a stand-in for the population. By resampling from your sample, you simulate what would happen if you could repeatedly sample from the population. It is remarkably effective even for small samples.

Hands-On: CLT with Allele Frequencies

Let us experience the Central Limit Theorem using realistic genetic data.

set_seed(42)
# Simulate allele frequency estimation from 1000 Genomes-style data

# True allele frequency of a common variant
let true_af = 0.23

# Simulate genotyping different numbers of individuals
let sample_sizes = [10, 30, 100, 500]

for n in sample_sizes {
    # Simulate 2000 studies, each genotyping n individuals
    let estimated_afs = []
    for study in 0..2000 {
        # Each individual contributes 2 alleles (diploid)
        let n_alleles = 2 * n
        let alt_count = rbinom(1, n_alleles, true_af) |> sum()
        let af_estimate = alt_count / n_alleles
        estimated_afs = estimated_afs + [af_estimate]
    }

    let se = stdev(estimated_afs)
    let theoretical_se = sqrt(true_af * (1.0 - true_af) / (2.0 * n))

    print("n={n}: SE={se:.4} (theoretical: {theoretical_se:.4})")
    histogram(estimated_afs, {bins: 40, title: "Allele Frequency Estimates (n={n})"})
}

# With n=10: estimates range wildly (0.05 to 0.50)
# With n=500: estimates tightly clustered around 0.23

This simulation shows exactly why GWAS studies need thousands of individuals. With 10 people, you cannot reliably estimate an allele frequency to better than ±10 percentage points. With 500, you can nail it to within ±1-2 percentage points.

Python and R Equivalents

Python:

import numpy as np
from scipy import stats

# Sampling distribution simulation
population = np.random.normal(125, 18, 100000)
sample_means = [np.mean(np.random.choice(population, 30)) for _ in range(1000)]
print(f"SE: {np.std(sample_means):.2f}")  # Should be close to 18/sqrt(30)

# Bootstrap
from scipy.stats import bootstrap
data = np.random.exponential(1, 50)
res = bootstrap((data,), np.median, n_resamples=5000)
print(f"95% CI: [{res.confidence_interval.low:.2f}, {res.confidence_interval.high:.2f}]")

# Standard error
se = np.std(data, ddof=1) / np.sqrt(len(data))

R:

# Sampling distribution
population <- rnorm(100000, mean = 125, sd = 18)
sample_means <- replicate(1000, mean(sample(population, 30)))
sd(sample_means)  # Empirical SE

# Bootstrap
library(boot)
boot_fn <- function(data, indices) median(data[indices])
results <- boot(data, boot_fn, R = 5000)
boot.ci(results, type = "perc")

# Standard error
se <- sd(data) / sqrt(length(data))

# Central Limit Theorem demo
par(mfrow = c(2, 2))
for (n in c(5, 10, 30, 100)) {
  means <- replicate(2000, mean(rexp(n, rate = 1)))
  hist(means, breaks = 40, main = paste("n =", n))
}

Exercises

Exercise 1: Experience the CLT

Take a uniform distribution (flat, definitely not normal) and show the CLT in action.

set_seed(42)
# Uniform population: values equally likely between 0 and 100
let uniform_pop = rnorm(100000, 50, 28.87)
# (Approximation — true uniform has SD = range/sqrt(12))

# TODO: Draw histograms of the population (should be flat-ish)
# TODO: Take 1000 samples of size n=5, compute means, draw histogram
# TODO: Repeat for n=10, n=30, n=100
# TODO: At what n does the sampling distribution look convincingly normal?
# TODO: Compute skewness at each n to quantify the convergence

Exercise 2: SE and Confidence

You measure tumor volumes in 25 mice (mean = 450 mm³, SD = 120 mm³).

let n = 25
let sample_mean = 450.0
let sample_sd = 120.0

# TODO: Compute the standard error
# TODO: Compute an approximate 95% CI using mean +/- 2*SE
# TODO: How large would n need to be for the 95% CI to have a width of +/-10 mm3?
# TODO: How large for +/-5 mm3?

Exercise 3: Bias Identification

For each scenario, identify the type of bias and explain how it could affect results.

1. A study of BRCA1 mutation frequency recruits subjects from a cancer genetics clinic.
2. A survival analysis of pancreatic cancer uses patients diagnosed 5+ years ago (all long-term survivors by definition).
3. RNA-seq libraries are prepared on two different days — all treatment samples on Day 1, all controls on Day 2.
4. A GWAS consortium publishes results for the 10 strongest associations and files the rest.

Exercise 4: Bootstrap a Correlation

Estimate the uncertainty in a correlation coefficient using the bootstrap.

set_seed(42)
# Gene expression vs. protein abundance (moderate correlation)
let n = 40
let gene_expr = rnorm(n, 5.0, 2.0)
let noise = rnorm(n, 0, 1.5)
let protein = gene_expr |> map(|x| 0.7 * x) |> zip(noise) |> map(|pair| pair.0 + pair.1)

let observed_r = cor(gene_expr, protein)
print("Observed correlation: {observed_r:.3}")

# TODO: Bootstrap the correlation 5000 times
# TODO: Compute the 95% bootstrap CI
# TODO: Is the correlation significantly different from zero?
# TODO: Plot the bootstrap distribution of r

Exercise 5: Power Simulation

Explore how effect size and variability interact with sample size.

set_seed(42)
# TODO: Run the trial simulation from the chapter, but now vary the effect size
# Test with differences of 2, 5, 10, and 20 mmHg (SD=18 throughout)
# At n=50 per arm, which effect sizes can you reliably detect?
# TODO: Now fix the difference at 5 mmHg and vary SD (10, 18, 30)
# How does variability affect the required sample size?

Key Takeaways

• Population vs. sample: You study a sample to learn about a population. The quality of inference depends on sample size and sampling method.
• The sampling distribution is the distribution of a statistic computed from repeated samples. It is narrower than the data distribution and centered at the true value.
• The Central Limit Theorem guarantees that sample means are approximately normal regardless of the population distribution, given sufficient sample size. This is why normal-based tests work so broadly.
• Standard error (SE = SD/√n) quantifies the precision of your estimate. Quadrupling n halves the SE.
• Bias (selection, survivorship, ascertainment, measurement, publication) is a systematic distortion that cannot be fixed by increasing n. Identify and prevent bias at the design stage.
• Statistical power is the probability of detecting a real effect. Underpowered studies waste resources and miss real effects. The four determinants of power are effect size, sample size, variability, and significance threshold.
• The bootstrap provides empirical estimates of standard errors and confidence intervals for any statistic, without relying on distributional assumptions.

What’s Next

You have now completed the foundations. You know how to summarize data (Day 2), understand its distributional shape (Day 3), reason about probabilities (Day 4), and appreciate the central role of sample size and sampling variability (Day 5). Starting next week, we put these foundations to work. Day 6 introduces confidence intervals — the formal framework for saying “I’m 95% sure the true value lies between here and here.” You will see how the standard error you learned today transforms into a rigorous statement about uncertainty, and why confidence intervals are more informative than p-values alone. The testing begins.

Day 6: Confidence Intervals — The Range of Truth

The Problem

Dr. Amara Chen’s pharmacology team has spent six months developing a novel kinase inhibitor for triple-negative breast cancer. After extensive optimization, they measure the half-maximal inhibitory concentration (IC50) across eight independent replicates: 11.2, 13.1, 12.8, 10.9, 14.2, 12.0, 11.7, and 12.5 nanomolar. The mean is 12.3 nM — an excellent result that would place their compound among the most potent in its class.

But when Dr. Chen presents these results to the medicinal chemistry team, the lead chemist asks the uncomfortable question: “If you ran the experiment again tomorrow, would you get 12.3 nM? Or could it be 15? Or 9?” The point estimate of 12.3 nM tells them where the center of their data is, but it says nothing about how confident they should be in that number. They need a range — a confidence interval — that captures the uncertainty inherent in measuring anything biological.

This chapter introduces the confidence interval: a range of plausible values for a population parameter, built from sample data. It is one of the most important and most misunderstood tools in all of biostatistics.

What Is a Confidence Interval?

Imagine you are trying to measure the height of a building, but your measuring tape is slightly stretchy. Each time you measure, you get a slightly different answer. A confidence interval is like saying: “Based on my eight measurements, I am 95% confident the true height is somewhere between 48.2 and 52.1 meters.”

More precisely: if you repeated your experiment 100 times and computed a 95% confidence interval each time, about 95 of those 100 intervals would contain the true population parameter. The remaining 5 would miss it entirely.

Common pitfall: A 95% CI does NOT mean “there is a 95% probability the true value is in this interval.” Once you compute a specific interval, the true value is either in it or it isn’t. The 95% refers to the procedure’s long-run success rate, not the probability for any single interval.

Point Estimates Are Not Enough

A point estimate is a single number — a sample mean, a proportion, a median. It is our best guess, but it carries no information about precision.

The confidence interval supplements the point estimate with a measure of uncertainty. Narrow intervals mean precise estimates; wide intervals mean the data leaves much room for doubt.

CI for a Mean: x-bar plus-or-minus t times SE

The most common confidence interval is for a population mean. The formula is:

CI = x-bar +/- t(alpha/2, df) x SE

Where:

• x-bar is the sample mean
• SE = s / sqrt(n) is the standard error of the mean
• t(alpha/2, df) is the critical value from the t-distribution with df = n - 1
• alpha = 1 - confidence level (for 95% CI, alpha = 0.05)

Why the t-Distribution for Small Samples?

When n is large (say, n > 30), the t-distribution closely resembles the normal distribution. But for small n — common in biology where each replicate is expensive — the t-distribution has heavier tails, producing wider intervals that honestly reflect our greater uncertainty.

Key insight: For biological experiments with n < 30, always use the t-distribution. Using z would give falsely narrow intervals that overstate your precision.

CI for a Proportion

When the variable is binary — mutation present/absent, responder/non-responder — we need a CI for a proportion p-hat = x/n.

Wald Interval (Simple but Flawed)

CI = p-hat +/- z x sqrt(p-hat(1 - p-hat) / n)

This is the textbook formula, but it performs poorly when p is near 0 or 1, or when n is small. It can even produce intervals that extend below 0 or above 1.

Wilson Interval (Preferred)

The Wilson score interval adjusts the center and width, and is recommended for most biological applications:

CI = (p-hat + z²/2n +/- z x sqrt(p-hat(1-p-hat)/n + z²/4n²)) / (1 + z²/n)

Clinical relevance: When reporting mutation carrier frequencies, drug response rates, or diagnostic sensitivity/specificity, always use Wilson intervals. Regulatory agencies expect intervals that behave properly even at extreme proportions.

CI for the Difference Between Two Means

Often the real question is not “what is the mean?” but “how much do two groups differ?” The CI for the difference between two independent means is:

CI = (x-bar1 - x-bar2) +/- t x SE_diff

Where SE_diff = sqrt(s1²/n1 + s2²/n2) for Welch’s approach.

The critical interpretation: If the CI for the difference includes zero, the data are consistent with no difference between the groups. If it excludes zero, the difference is statistically significant.

Bootstrap Confidence Intervals

What if your statistic is a median, a ratio, or something with no tidy formula? The bootstrap is a computer-intensive method that works for any statistic:

1. Resample your data with replacement, same size as original
2. Compute the statistic on the resample
3. Repeat 10,000 times
4. Take the 2.5th and 97.5th percentiles of the bootstrap distribution

This is called the percentile method. No assumptions about normality or distribution shape are required.

Key insight: Bootstrap CIs are the Swiss army knife of interval estimation. When in doubt, bootstrap it.

What Controls CI Width?

Three factors determine how wide or narrow your confidence interval will be:

This is why power calculations matter: before an experiment, you choose n to achieve a CI narrow enough to be scientifically useful.

Confidence Intervals in BioLang

IC50 Confidence Interval — Parametric

# IC50 measurements (nM) from 8 replicates
let ic50 = [11.2, 13.1, 12.8, 10.9, 14.2, 12.0, 11.7, 12.5]

let n = len(ic50)
let x_bar = mean(ic50)
let se = stdev(ic50) / sqrt(n)

# 95% CI using normal approximation (for small n, t > z)
let t_crit = qnorm(0.975)
let ci_lower = x_bar - t_crit * se
let ci_upper = x_bar + t_crit * se

print("IC50 mean: {x_bar:.2} nM")
print("95% CI: [{ci_lower:.2}, {ci_upper:.2}] nM")
print("Standard error: {se:.3} nM")
print("Critical value: {t_crit:.3}")

IC50 Confidence Interval — Bootstrap

set_seed(42)
# Bootstrap CI: no distributional assumptions

let ic50 = [11.2, 13.1, 12.8, 10.9, 14.2, 12.0, 11.7, 12.5]

# Bootstrap: resample 10,000 times, compute mean each time
let n_boot = 10000
let boot_means = []
for i in range(0, n_boot) {
    let resample = []
    for j in range(0, len(ic50)) {
        resample = append(resample, ic50[random_int(0, len(ic50) - 1)])
    }
    boot_means = append(boot_means, mean(resample))
}

# Percentile method
let boot_lower = quantile(boot_means, 0.025)
let boot_upper = quantile(boot_means, 0.975)

print("Bootstrap 95% CI: [{boot_lower:.2}, {boot_upper:.2}] nM")

# Visualize the bootstrap distribution
histogram(boot_means, {bins: 50, title: "Bootstrap Distribution of IC50 Mean", x_label: "Mean IC50 (nM)"})

Bootstrap CI for Median (No Parametric Formula Exists)

set_seed(42)
# Gene expression values (FPKM) — skewed distribution
let expression = [0.1, 0.3, 0.8, 1.2, 1.5, 2.1, 3.4, 8.7, 12.1, 45.6]

let obs_median = median(expression)

# Bootstrap the median
let n_boot = 10000
let boot_medians = []
for i in range(0, n_boot) {
    let resample = []
    for j in range(0, len(expression)) {
        resample = append(resample, expression[random_int(0, len(expression) - 1)])
    }
    boot_medians = append(boot_medians, median(resample))
}
let ci_lower = quantile(boot_medians, 0.025)
let ci_upper = quantile(boot_medians, 0.975)

print("Observed median: {obs_median:.2} FPKM")
print("Bootstrap 95% CI for median: [{ci_lower:.2}, {ci_upper:.2}] FPKM")

Error Bar Plot: Comparing Drug Concentrations

# IC50 values for three drug candidates
let drug_a = [12.3, 11.8, 13.1, 12.0, 11.5, 12.7, 13.4, 12.1]
let drug_b = [25.1, 28.3, 22.7, 26.9, 24.5, 27.1, 23.8, 25.6]
let drug_c = [8.2, 9.1, 7.5, 8.8, 10.2, 8.0, 9.5, 7.8]

let drugs = ["Drug A", "Drug B", "Drug C"]
let means = [mean(drug_a), mean(drug_b), mean(drug_c)]

# Compute 95% CIs for each
let compute_ci = |data| {
  let n = len(data)
  let se = stdev(data) / sqrt(n)
  let t_crit = qnorm(0.975)
  [mean(data) - t_crit * se, mean(data) + t_crit * se]
}

let ci_a = compute_ci(drug_a)
let ci_b = compute_ci(drug_b)
let ci_c = compute_ci(drug_c)

print("Drug A: {means[0]:.1} nM, 95% CI [{ci_a[0]:.1}, {ci_a[1]:.1}]")
print("Drug B: {means[1]:.1} nM, 95% CI [{ci_b[0]:.1}, {ci_b[1]:.1}]")
print("Drug C: {means[2]:.1} nM, 95% CI [{ci_c[0]:.1}, {ci_c[1]:.1}]")

# Bar chart with error bars
bar_chart(drugs, means, {title: "IC50 Comparison with 95% CIs", y_label: "IC50 (nM)", error_bars: [ci_a, ci_b, ci_c]})

CI for Difference Between Two Means

# Compare tumor volume between treated and control mice
let treated = [180, 210, 165, 225, 195, 172, 218, 198]
let control = [485, 512, 468, 530, 495, 478, 521, 503]

let diff = mean(treated) - mean(control)
let se_diff = sqrt(variance(treated) / len(treated) + variance(control) / len(control))
let df = len(treated) + len(control) - 2
let t_crit = qnorm(0.975)  # approximate for moderate df

let ci_lower = diff - t_crit * se_diff
let ci_upper = diff + t_crit * se_diff

print("Mean difference: {diff:.1} mm^3")
print("95% CI for difference: [{ci_lower:.1}, {ci_upper:.1}] mm^3")

if ci_upper < 0 {
  print("CI excludes zero: treatment significantly reduces tumor volume")
} else {
  print("CI includes zero: cannot rule out no difference")
}

Vaccine Efficacy CI (Proportion)

set_seed(42)
# Clinical trial: 15 of 200 vaccinated got infected vs 60 of 200 placebo
let p_vacc = 15 / 200
let p_plac = 60 / 200
let efficacy = 1.0 - (p_vacc / p_plac)

print("Vaccine efficacy: {efficacy * 100:.1}%")

# CI for proportion (vaccinated group infection rate)
let n = 200
let z = 1.96
let se_p = sqrt(p_vacc * (1.0 - p_vacc) / n)
let ci_lower_p = p_vacc - z * se_p
let ci_upper_p = p_vacc + z * se_p

print("Infection rate (vaccinated): {p_vacc*100:.1}%")
print("95% CI for infection rate: [{ci_lower_p*100:.1}%, {ci_upper_p*100:.1}%]")

# Bootstrap CI for vaccine efficacy itself
let vacc_outcomes = flatten([repeat(1, 15), repeat(0, 185)])
let plac_outcomes = flatten([repeat(1, 60), repeat(0, 140)])

let n_boot = 10000
let boot_eff = []
for i in range(0, n_boot) {
    let v_resample = []
    let p_resample = []
    for j in range(0, len(vacc_outcomes)) {
        v_resample = append(v_resample, vacc_outcomes[random_int(0, len(vacc_outcomes) - 1)])
        p_resample = append(p_resample, plac_outcomes[random_int(0, len(plac_outcomes) - 1)])
    }
    let pv = mean(v_resample)
    let pp = mean(p_resample)
    let eff = if pp == 0.0 then 0.0 else 1.0 - (pv / pp)
    boot_eff = append(boot_eff, eff)
}

let eff_ci = [quantile(boot_eff, 0.025), quantile(boot_eff, 0.975)]
print("Bootstrap 95% CI for efficacy: [{eff_ci[0]*100:.1}%, {eff_ci[1]*100:.1}%]")

Python:

import numpy as np
from scipy import stats

ic50 = [11.2, 13.1, 12.8, 10.9, 14.2, 12.0, 11.7, 12.5]
ci = stats.t.interval(0.95, df=len(ic50)-1,
                       loc=np.mean(ic50),
                       scale=stats.sem(ic50))
print(f"95% CI: [{ci[0]:.2f}, {ci[1]:.2f}]")

# Bootstrap
boot = [np.mean(np.random.choice(ic50, len(ic50))) for _ in range(10000)]
print(f"Bootstrap CI: [{np.percentile(boot, 2.5):.2f}, {np.percentile(boot, 97.5):.2f}]")

R:

ic50 <- c(11.2, 13.1, 12.8, 10.9, 14.2, 12.0, 11.7, 12.5)
t.test(ic50)$conf.int

# Bootstrap
library(boot)
boot_fn <- function(data, i) mean(data[i])
b <- boot(ic50, boot_fn, R = 10000)
boot.ci(b, type = "perc")

Exercises

Exercise 1: Compute a CI by Hand

Eight mice on a high-fat diet had cholesterol levels: 215, 228, 197, 241, 209, 233, 220, 212 mg/dL. Compute the 95% CI for the mean cholesterol.

let cholesterol = [215, 228, 197, 241, 209, 233, 220, 212]

# TODO: Compute mean, SE, critical value, and the 95% CI
# Hint: df = n - 1, use qnorm(0.975) as approximate critical value

Exercise 2: Bootstrap a Ratio

Gene A has FPKM values [2.1, 3.4, 1.8, 4.2, 2.9] in tumor and [1.0, 1.2, 0.9, 1.5, 1.1] in normal. Bootstrap a 95% CI for the tumor/normal fold change of medians.

let tumor = [2.1, 3.4, 1.8, 4.2, 2.9]
let normal = [1.0, 1.2, 0.9, 1.5, 1.1]

# TODO: Bootstrap the ratio median(tumor) / median(normal)
# Use n_boot = 10000, then extract 2.5th and 97.5th percentiles with quantile()

Exercise 3: Overlapping CIs

Compute 95% CIs for Drug X (IC50: [5.2, 6.1, 4.8, 5.5, 6.3, 5.0]) and Drug Y (IC50: [5.8, 6.5, 7.2, 6.0, 5.9, 6.8]). Do the CIs overlap? What does this suggest?

let drug_x = [5.2, 6.1, 4.8, 5.5, 6.3, 5.0]
let drug_y = [5.8, 6.5, 7.2, 6.0, 5.9, 6.8]

# TODO: Compute CIs for both, then compute CI for the difference
# Note: overlapping CIs do NOT necessarily mean non-significant difference

Exercise 4: Effect of Sample Size

Starting with n = 5 replicates drawn from the IC50 data, increase to n = 10, 20, 50, and 100 (use bootstrap resampling to simulate larger samples). Plot CI width vs sample size.

let ic50 = [11.2, 13.1, 12.8, 10.9, 14.2, 12.0, 11.7, 12.5]

# TODO: For each sample size, bootstrap to simulate, compute CI width
# Plot sample size vs CI width using line_plot

Key Takeaways

• A confidence interval gives a range of plausible values for a population parameter, not just a point estimate
• The 95% in “95% CI” refers to the long-run coverage rate of the procedure, not the probability for a specific interval
• For small samples (n < 30), always use the t-distribution — it accounts for extra uncertainty
• Bootstrap CIs work for any statistic (median, ratio, fold change) without distributional assumptions
• CI width shrinks with larger n, lower variability, and lower confidence level
• A CI for the difference that includes zero means the data are consistent with no difference
• CIs are more informative than p-values alone: they tell you both significance AND the plausible magnitude of an effect

What’s Next

Tomorrow we formalize the logic behind “ruling out chance” with hypothesis testing. You will learn to frame biological questions as null and alternative hypotheses, compute p-values, and understand the courtroom analogy that makes the whole framework click. Confidence intervals and hypothesis tests are two sides of the same coin — a 95% CI that excludes zero corresponds exactly to a p-value less than 0.05.

Day 7: Hypothesis Testing — Asking Precise Questions

The Problem

Dr. Kenji Nakamura has spent three years developing a blood-based biomarker panel for early Alzheimer’s detection. His team measures plasma levels of phosphorylated tau (p-tau217) in 40 cognitively normal individuals and 40 patients with confirmed early-stage Alzheimer’s. The mean p-tau217 level in the Alzheimer’s group is 3.8 pg/mL, compared to 2.9 pg/mL in controls. The difference looks promising — nearly 30% higher.

But when Dr. Nakamura submits to the FDA for breakthrough device designation, the reviewer’s response is blunt: “Your biomarker shows a numerical difference. Can you demonstrate this isn’t just sampling noise? What is the probability of seeing a difference this large if the biomarker has no real diagnostic value?” This is the fundamental question that hypothesis testing answers.

The stakes are enormous. If the biomarker works, millions of patients could be diagnosed years earlier, when interventions are most effective. If it doesn’t — if the observed difference is just statistical noise — pursuing it wastes hundreds of millions in development costs and, worse, could lead to false diagnoses.

What Is Hypothesis Testing?

Think of hypothesis testing as a courtroom trial for your scientific claim.

• The defendant is the null hypothesis (H0): “There is no effect.” In the courtroom, the defendant is presumed innocent.
• The prosecution’s evidence is your data. You are trying to show the evidence is so overwhelming that the “innocence” explanation is implausible.
• The verdict is either “guilty” (reject H0) or “not proven” (fail to reject H0). Notice: the jury never declares the defendant “innocent” — just that the evidence was insufficient.

Key insight: Hypothesis testing never proves your theory is true. It only tells you whether the data are inconsistent enough with “no effect” that you can reject that explanation with a specified level of confidence.

The Five Steps of Hypothesis Testing

The Null Hypothesis (H0)

The null hypothesis is the “boring” explanation — the default assumption of no effect, no difference, no relationship. It is what you assume until the data force you to abandon it.

The Alternative Hypothesis (H1)

The alternative is what you actually believe — the “interesting” claim.

• Two-tailed: H1: mu1 != mu2 (the groups differ in either direction)
• One-tailed: H1: mu1 > mu2 (specifically higher) or H1: mu1 < mu2 (specifically lower)

Common pitfall: Do not choose one-tailed vs two-tailed after looking at your data. This decision must be made before the experiment, based on your scientific question. Switching from two-tailed to one-tailed after seeing results halves your p-value — that is scientific fraud.

The p-Value: Most Misunderstood Number in Science

The p-value is the probability of observing data as extreme as (or more extreme than) what you got, assuming H0 is true.

What the p-value IS:

• A measure of how surprising your data are under the null hypothesis
• A continuous measure of evidence — smaller p = more evidence against H0
• The probability of the data given H0: P(data | H0)

What the p-value IS NOT:

• The probability that H0 is true: NOT P(H0 | data)
• The probability your result is due to chance
• The probability of making an error
• A measure of effect size or practical importance

Type I and Type II Errors

Every decision carries the risk of being wrong:

• Type I error (alpha): You claim a drug works when it doesn’t. A patient receives ineffective treatment.
• Type II error (beta): You miss a real effect. An effective drug gets shelved.

Clinical relevance: In drug safety testing, alpha is typically set very low (0.01 or even 0.001) because a Type I error means approving a dangerous drug. In exploratory genomics, higher alpha (0.05 or even 0.10) is acceptable because you will validate hits in follow-up experiments.

Statistical vs Practical Significance

A p-value of 0.001 does not mean the effect is large or important. With enough data, trivially small effects become “statistically significant.”

Always report effect sizes alongside p-values. We will dedicate Day 19 entirely to this topic.

The z-Test: The Simplest Hypothesis Test

When the population standard deviation sigma is known (rare, but a good starting point), the z-test compares a sample mean to a known value:

z = (x-bar - mu0) / (sigma / sqrt(n))

Under H0, z follows a standard normal distribution N(0, 1).

One-Tailed vs Two-Tailed Tests

Hypothesis Testing in BioLang

z-Test on Biomarker Data

# Alzheimer's biomarker study
# Known population SD from large reference database: sigma = 1.2 pg/mL
# Expected normal level: mu0 = 2.9 pg/mL
# Observed in 40 Alzheimer's patients:
let ad_levels = [3.2, 4.1, 3.8, 2.7, 4.5, 3.3, 3.9, 4.2, 3.1, 3.6,
                 4.0, 3.5, 4.3, 3.7, 2.9, 3.8, 4.1, 3.4, 3.6, 4.4,
                 3.0, 3.9, 4.2, 3.3, 3.7, 4.0, 3.5, 3.8, 4.1, 3.2,
                 3.6, 3.4, 4.3, 3.1, 3.9, 3.7, 4.0, 3.5, 3.8, 3.3]

# Compute z-test manually (known sigma)
let n = len(ad_levels)
let x_bar = mean(ad_levels)
let se = 1.2 / sqrt(n)
let z_stat = (x_bar - 2.9) / se
let p_val = 2.0 * (1.0 - pnorm(abs(z_stat), 0, 1))

print("z-statistic: {z_stat:.4}")
print("p-value (two-tailed): {p_val:.6}")
print("Mean observed: {x_bar:.3} pg/mL")

if p_val < 0.05 {
  print("Reject H0: Alzheimer's group significantly differs from normal level")
} else {
  print("Fail to reject H0: Insufficient evidence of difference")
}

Visualizing the Null Distribution

# Show where our test statistic falls on the null distribution
let z_obs = 4.71  # from the z-test above

# Generate the null distribution (standard normal)
let x_vals = range(-4.0, 4.0, 0.01)
let y_vals = x_vals |> map(|x| dnorm(x, 0, 1))

# Plot the null distribution with our observed z marked
density(x_vals, {title: "Null Distribution (Standard Normal)", x_label: "z-statistic", y_label: "Density", vlines: [z_obs], shade_above: 1.96, shade_below: -1.96})

print("Critical value (two-tailed, alpha=0.05): +/- 1.96")
print("Our z = {z_obs} falls far in the rejection region")

Binomial Test: Is This Mutation Rate Elevated?

# In a cancer cohort, 18 of 100 patients carry a specific BRCA2 variant
# Population frequency is known to be 8%
# Compute binomial test using dbinom: P(X >= 18) when X ~ Binom(100, 0.08)
let p_val = 0.0
for k in range(18, 101) {
    p_val = p_val + dbinom(k, 100, 0.08)
}

print("Observed proportion: 18/100 = 18%")
print("Expected under H0: 8%")
print("p-value (one-sided): {p_val:.6}")

if p_val < 0.05 {
  print("Reject H0: Mutation rate is significantly elevated in this cohort")
} else {
  print("Fail to reject H0")
}

Complete Hypothesis Test Workflow

# Full workflow: Is mean platelet count elevated in a disease cohort?
# Reference population: mu = 250 (x10^3/uL), sigma = 50
# Our 25 patients:
let platelets = [280, 310, 265, 295, 275, 320, 290, 305, 260, 285,
                 300, 270, 315, 288, 292, 278, 308, 282, 298, 272,
                 310, 295, 268, 302, 288]

# Step 1: State hypotheses
print("H0: mu = 250 (platelet count is normal)")
print("H1: mu > 250 (platelet count is elevated)")
print("alpha = 0.05, one-tailed test\n")

# Step 2: Compute test statistic
let n = len(platelets)
let x_bar = mean(platelets)
let se = 50 / sqrt(n)  # sigma is known
let z = (x_bar - 250) / se

# Step 3: Find p-value (one-tailed)
let p = 1.0 - pnorm(z, 0, 1)

# Step 4: Decision
print("Sample mean: {x_bar:.1}")
print("z-statistic: {z:.4}")
print("p-value (one-tailed): {p:.6}")

if p < 0.05 {
  print("\nDecision: Reject H0 at alpha = 0.05")
  print("Conclusion: Platelet count is significantly elevated in this cohort")
} else {
  print("\nDecision: Fail to reject H0")
}

Interpreting p-Values with Simulated Data

set_seed(42)
# Demonstrate: under H0 (no effect), p-values are uniformly distributed

let p_values = []
for i in 1..1000 {
  # Generate two samples from the SAME distribution (H0 is true)
  let group1 = rnorm(20, 10, 3)
  let group2 = rnorm(20, 10, 3)
  let z_stat = (mean(group1) - mean(group2)) / (3.0 / sqrt(20))
  let p_val = 2.0 * (1.0 - pnorm(abs(z_stat), 0, 1))
  p_values = append(p_values, p_val)
}

# Count false positives at alpha = 0.05
let false_pos = p_values |> filter(|p| p < 0.05) |> len()
print("False positives out of 1000 null tests: {false_pos}")
print("Expected: ~50 (5% of 1000)")

histogram(p_values, {title: "p-Value Distribution Under the Null", x_label: "p-value", bins: 20})

Connecting CIs and Hypothesis Tests

# Demonstrate: a 95% CI that excludes the null value corresponds to p < 0.05
let ad_levels = [3.2, 4.1, 3.8, 2.7, 4.5, 3.3, 3.9, 4.2, 3.1, 3.6,
                 4.0, 3.5, 4.3, 3.7, 2.9, 3.8, 4.1, 3.4, 3.6, 4.4]

# z-test: is the mean different from 2.9?
let n = len(ad_levels)
let x_bar = mean(ad_levels)
let se = 1.2 / sqrt(n)
let z_stat = (x_bar - 2.9) / se
let z_p = 2.0 * (1.0 - pnorm(abs(z_stat), 0, 1))
print("z-test p-value: {z_p:.6}")

# 95% CI for the mean (using known sigma)
let se2 = 1.2 / sqrt(n)
let ci_lower = x_bar - 1.96 * se2
let ci_upper = x_bar + 1.96 * se2
print("95% CI: [{ci_lower:.3}, {ci_upper:.3}]")
print("Null value (2.9) is {'outside' if ci_lower > 2.9 or ci_upper < 2.9 else 'inside'} the CI")
print("This matches the hypothesis test: p < 0.05 <=> CI excludes null value")

Python:

import numpy as np
from scipy import stats

ad_levels = [3.2, 4.1, 3.8, 2.7, 4.5, 3.3, 3.9, 4.2, 3.1, 3.6,
             4.0, 3.5, 4.3, 3.7, 2.9, 3.8, 4.1, 3.4, 3.6, 4.4,
             3.0, 3.9, 4.2, 3.3, 3.7, 4.0, 3.5, 3.8, 4.1, 3.2,
             3.6, 3.4, 4.3, 3.1, 3.9, 3.7, 4.0, 3.5, 3.8, 3.3]

# z-test (manual — scipy doesn't have a built-in z-test for means)
z = (np.mean(ad_levels) - 2.9) / (1.2 / np.sqrt(len(ad_levels)))
p = 2 * (1 - stats.norm.cdf(abs(z)))
print(f"z = {z:.4f}, p = {p:.6f}")

# Binomial test
result = stats.binomtest(18, 100, 0.08, alternative='greater')
print(f"Binomial test p = {result.pvalue:.6f}")

R:

ad_levels <- c(3.2, 4.1, 3.8, 2.7, 4.5, 3.3, 3.9, 4.2, 3.1, 3.6,
               4.0, 3.5, 4.3, 3.7, 2.9, 3.8, 4.1, 3.4, 3.6, 4.4,
               3.0, 3.9, 4.2, 3.3, 3.7, 4.0, 3.5, 3.8, 4.1, 3.2,
               3.6, 3.4, 4.3, 3.1, 3.9, 3.7, 4.0, 3.5, 3.8, 3.3)

# z-test (using BSDA package, or manual)
z <- (mean(ad_levels) - 2.9) / (1.2 / sqrt(length(ad_levels)))
p <- 2 * pnorm(-abs(z))
cat(sprintf("z = %.4f, p = %.6f\n", z, p))

# Binomial test
binom.test(18, 100, p = 0.08, alternative = "greater")

Exercises

Exercise 1: Formulate Hypotheses

For each scenario, write the null and alternative hypotheses. State whether you would use a one-tailed or two-tailed test and why.

a) Does a new antibiotic reduce bacterial colony counts compared to placebo? b) Is the GC content of a newly sequenced genome different from the expected 42%? c) Do patients with the variant allele have higher LDL cholesterol?

Exercise 2: z-Test on Gene Expression

A reference database reports the mean expression of housekeeping gene GAPDH as 8.5 log2-CPM with sigma = 0.8 across thousands of samples. Your RNA-seq experiment on 15 samples yields a mean of 7.9. Is your experiment’s GAPDH level significantly different?

let gapdh_expression = [7.5, 8.1, 7.8, 7.6, 8.3, 7.2, 8.0, 7.9,
                        8.2, 7.4, 7.7, 8.1, 7.6, 8.4, 7.3]

# TODO: Perform z-test with mu=8.5, sigma=0.8
# TODO: Interpret the result — what might explain a difference?

Exercise 3: Simulate Type I Error Rate

Run 10,000 z-tests where H0 is true (both groups from the same distribution). Count what fraction of p-values fall below 0.05, 0.01, and 0.001. Do these match expectations?

# TODO: Simulate 10,000 null tests
# TODO: Count p < 0.05, p < 0.01, p < 0.001
# TODO: Compare to expected rates (5%, 1%, 0.1%)

Exercise 4: One-Tailed vs Two-Tailed

Using the Alzheimer’s biomarker data, compute the p-value for both a one-tailed test (H1: AD levels are higher) and a two-tailed test. What is the relationship between the two p-values?

let ad_levels = [3.2, 4.1, 3.8, 2.7, 4.5, 3.3, 3.9, 4.2, 3.1, 3.6,
                 4.0, 3.5, 4.3, 3.7, 2.9, 3.8, 4.1, 3.4, 3.6, 4.4]

# TODO: Compute z-stat manually, then get two-tailed and one-tailed p-values
# Two-tailed: 2 * (1 - pnorm(abs(z), 0, 1))
# One-tailed: 1 - pnorm(z, 0, 1)
# TODO: What is the mathematical relationship?

Key Takeaways

• Hypothesis testing uses the courtroom analogy: H0 (innocence) is assumed until the evidence (data) is overwhelming
• The p-value is the probability of data this extreme under H0 — it is NOT the probability H0 is true
• Type I error (false positive) is controlled by alpha; Type II error (false negative) is controlled by power
• Statistical significance (small p) does not imply practical significance (large effect)
• The z-test is the simplest hypothesis test, applicable when sigma is known
• Always state hypotheses and choose alpha before looking at data
• Under the null, p-values are uniformly distributed: at alpha = 0.05, exactly 5% of null tests will be “significant” by chance

What’s Next

Tomorrow we move from the z-test (which requires known sigma) to the workhorse of biological research: the t-test. You will learn independent, paired, and Welch’s versions, check assumptions with Shapiro-Wilk and Levene’s tests, and quantify effect sizes with Cohen’s d. If hypothesis testing is the question, the t-test is the answer for two-group comparisons.

Day 8: Comparing Two Groups — The t-Test

The Problem

Dr. Sofia Reyes is a cancer biologist studying BRCA1 expression in breast tissue. She has RNA-seq data from 12 tumor samples and 12 matched normal samples from the same patients. The mean BRCA1 expression in tumors is 4.2 log2-CPM versus 6.8 log2-CPM in normals — a 2.6-fold reduction. But with only 12 samples per group and considerable biological variability, can she confidently claim BRCA1 is downregulated in tumors?

She cannot use a z-test because the population standard deviation is unknown — she must estimate it from the data itself. She needs the t-test, the most widely used statistical test in biomedical research. But which version? Her samples are paired (tumor and normal from the same patient), which adds another consideration. And before running any test, she should verify that the data meet the test’s assumptions.

This chapter covers the t-test in all its forms: independent, Welch’s, paired, and one-sample. You will learn when each is appropriate, how to check assumptions, and how to quantify the magnitude of differences with Cohen’s d.

What Is the t-Test?

The t-test asks: “Is the difference between two group means larger than what we would expect from random sampling variation alone?”

Think of it this way: you have two piles of measurements. The t-test weighs how far apart the piles’ centers are, relative to how spread out each pile is. If the piles are far apart and tight, the difference is convincing. If they overlap substantially, it is not.

The t-statistic = (difference in means) / (standard error of the difference)

A larger t-statistic means more evidence of a real difference.

The Four Flavors of t-Test

Independent Two-Sample t-Test

Assumptions

1. Independence: Observations within and between groups are independent
2. Normality: Data in each group are approximately normally distributed
3. Equal variances: Both groups have similar spread (homoscedasticity)

The Pooled Standard Error

When variances are assumed equal, we pool them for a better estimate:

s_p = sqrt(((n1-1)s1^2 + (n2-1)s2^2) / (n1 + n2 - 2))

Degrees of freedom: df = n1 + n2 - 2

Welch’s t-Test: The Safer Default

Welch’s t-test does not assume equal variances. It uses each group’s own variance estimate and adjusts the degrees of freedom downward with the Welch-Satterthwaite equation.

Key insight: Welch’s t-test is almost always the better default choice. It performs nearly as well as the pooled t-test when variances ARE equal, and much better when they are not. Most modern statistical software (including R’s t.test()) uses Welch’s version by default.

Paired t-Test: Matched Samples

When observations are naturally paired — tumor/normal from the same patient, before/after treatment on the same subject — the paired t-test is far more powerful because it controls for inter-subject variability.

The trick: compute the difference for each pair, then perform a one-sample t-test on the differences:

t = d-bar / (s_d / sqrt(n))

Where d-bar is the mean of the paired differences and s_d is their standard deviation.

Common pitfall: Using an independent t-test on paired data wastes statistical power. If you have natural pairs, always use the paired test. Conversely, using a paired test on unpaired data gives wrong results.

Checking Assumptions

Normality: Shapiro-Wilk Test

The Shapiro-Wilk test checks whether data could have come from a normal distribution.

• H0: Data are normally distributed
• If p > 0.05, normality assumption is reasonable
• If p < 0.05, data are significantly non-normal

Also use QQ plots: if points fall along the diagonal line, data are approximately normal.

Equal Variances: Levene’s Test

Levene’s test checks whether two groups have equal variances.

• H0: Variances are equal
• If p > 0.05, equal variance assumption is reasonable
• If p < 0.05, use Welch’s t-test (or just always use Welch’s)

Cohen’s d: Quantifying Effect Size

A p-value tells you whether a difference exists. Cohen’s d tells you how large it is, in standard deviation units:

d = (x-bar1 - x-bar2) / s_pooled

Key insight: A large p-value with a large Cohen’s d suggests you are underpowered — you may have a real effect but too few samples to detect it. A small p-value with a tiny Cohen’s d suggests the effect, while real, may not be biologically meaningful.

The t-Test in BioLang

Independent Two-Sample t-Test: Gene Expression

# BRCA1 expression (log2-CPM) in tumor vs normal breast tissue
let tumor  = [3.8, 4.5, 4.1, 3.2, 4.8, 3.9, 4.3, 5.1, 3.6, 4.0, 4.7, 3.5]
let normal = [6.2, 7.1, 6.5, 7.4, 6.8, 7.0, 6.3, 7.2, 6.9, 6.6, 7.3, 6.1]

# Default: Welch's t-test (unequal variances)
let result = ttest(tumor, normal)
print("=== Welch's t-test: BRCA1 Tumor vs Normal ===")
print("t-statistic: {result.statistic:.4}")
print("p-value: {result.p_value:.2e}")
print("Degrees of freedom: {result.df:.1}")
print("Mean tumor: {mean(tumor):.2}, Mean normal: {mean(normal):.2}")
print("Difference: {mean(tumor) - mean(normal):.2} log2-CPM")

# Effect size (Cohen's d inline)
let d = (mean(tumor) - mean(normal)) / sqrt((variance(tumor) + variance(normal)) / 2.0)
print("Cohen's d: {d:.3}")

# Visualize
let bp_table = table({"Tumor": tumor, "Normal": normal})
boxplot(bp_table, {title: "BRCA1 Expression: Tumor vs Normal"})

Checking Assumptions

let tumor  = [3.8, 4.5, 4.1, 3.2, 4.8, 3.9, 4.3, 5.1, 3.6, 4.0, 4.7, 3.5]
let normal = [6.2, 7.1, 6.5, 7.4, 6.8, 7.0, 6.3, 7.2, 6.9, 6.6, 7.3, 6.1]

# 1. Normality check: use QQ plots for visual assessment
# (no built-in Shapiro-Wilk; use QQ plots + summary stats)
let s_tumor = summary(tumor)
let s_normal = summary(normal)
print("Tumor summary:  {s_tumor}")
print("Normal summary: {s_normal}")

# 2. Equal variance check: compare variances from summary()
let var_ratio = variance(tumor) / variance(normal)
print("Variance ratio (tumor/normal): {var_ratio:.3}")
if var_ratio > 2.0 or var_ratio < 0.5 {
  print("Variances appear unequal -> use Welch's t-test (the default)")
} else {
  print("Variances appear similar -> pooled t-test is also valid")
}

# 3. QQ plots for visual normality assessment
qq_plot(tumor, {title: "QQ Plot: Tumor BRCA1 Expression"})
qq_plot(normal, {title: "QQ Plot: Normal BRCA1 Expression"})

Paired t-Test: Before/After Treatment

# Tumor volume (mm^3) before and after 6 weeks of treatment
# Same 10 patients measured at both time points
let before = [245, 312, 198, 367, 289, 421, 156, 334, 278, 305]
let after  = [180, 245, 165, 298, 220, 350, 132, 270, 210, 248]

# Paired t-test: accounts for patient-to-patient variability
let result = ttest_paired(before, after)
print("=== Paired t-test: Tumor Volume Before vs After Treatment ===")
print("t-statistic: {result.statistic:.4}")
print("p-value: {result.p_value:.6}")

# Show the paired differences
let diffs = zip(before, after) |> map(|pair| pair[0] - pair[1])
print("Mean reduction: {mean(diffs):.1} mm^3")
print("Individual reductions: {diffs}")

# Compare: what if we wrongly used an independent t-test?
let wrong_result = ttest(before, after)
print("\nWrong (independent) t-test p-value: {wrong_result.p_value:.6}")
print("Correct (paired) t-test p-value: {result.p_value:.6}")
print("Paired test is more powerful because it removes inter-patient variability")

# Visualize paired differences
histogram(diffs, {title: "Distribution of Tumor Volume Reductions", x_label: "Reduction (mm^3)", bins: 8})

One-Sample t-Test

# Is the GC content of our assembled genome different from the expected 41%?
let gc_per_contig = [40.2, 41.5, 39.8, 42.1, 40.7, 41.3, 39.5, 42.4,
                     40.1, 41.8, 40.5, 41.0, 39.9, 41.6, 40.3]

let result = ttest_one(gc_per_contig, 41.0)
print("=== One-sample t-test: GC Content vs Expected 41% ===")
print("Sample mean: {mean(gc_per_contig):.2}%")
print("t-statistic: {result.statistic:.4}")
print("p-value: {result.p_value:.4}")

if result.p_value > 0.05 {
  print("No significant deviation from expected GC content")
}

Complete Workflow: Multiple Genes

# Test multiple genes at once
let genes = ["BRCA1", "TP53", "MYC", "GAPDH", "EGFR"]

let tumor_expr = [
  [3.8, 4.5, 4.1, 3.2, 4.8, 3.9, 4.3, 5.1, 3.6, 4.0, 4.7, 3.5],
  [2.1, 1.8, 2.5, 1.4, 2.2, 1.9, 2.0, 1.6, 2.3, 1.7, 2.4, 1.5],
  [9.2, 10.1, 8.8, 9.5, 10.3, 9.7, 8.6, 9.9, 10.5, 9.1, 9.8, 10.2],
  [8.1, 8.3, 7.9, 8.2, 8.0, 8.4, 7.8, 8.1, 8.3, 8.0, 8.2, 7.9],
  [7.5, 8.2, 7.8, 8.5, 7.1, 8.0, 7.6, 8.3, 7.9, 8.1, 7.4, 8.4]
]

let normal_expr = [
  [6.2, 7.1, 6.5, 7.4, 6.8, 7.0, 6.3, 7.2, 6.9, 6.6, 7.3, 6.1],
  [5.8, 6.2, 5.5, 6.0, 5.9, 6.3, 5.7, 6.1, 5.6, 6.4, 5.8, 6.0],
  [5.1, 5.4, 4.9, 5.3, 5.6, 5.0, 5.2, 5.5, 4.8, 5.7, 5.1, 5.4],
  [8.0, 8.2, 7.8, 8.3, 8.1, 8.0, 8.4, 7.9, 8.2, 8.1, 8.3, 8.0],
  [5.0, 5.3, 4.8, 5.1, 5.5, 4.9, 5.2, 5.4, 4.7, 5.6, 5.0, 5.3]
]

print("Gene       | t-stat | p-value    | Cohen's d | Interpretation")
print("-----------|--------|------------|-----------|---------------")

for i in 0..len(genes) {
  let result = ttest(tumor_expr[i], normal_expr[i])
  let d = (mean(tumor_expr[i]) - mean(normal_expr[i])) / sqrt((variance(tumor_expr[i]) + variance(normal_expr[i])) / 2.0)
  let interp = if abs(d) > 0.8 then "Large" else if abs(d) > 0.5 then "Medium" else "Small"
  print("{genes[i]:<10} | {result.statistic:>6.2} | {result.p_value:>10.2e} | {d:>9.3} | {interp}")
}

Python:

from scipy import stats
import numpy as np

tumor  = [3.8, 4.5, 4.1, 3.2, 4.8, 3.9, 4.3, 5.1, 3.6, 4.0, 4.7, 3.5]
normal = [6.2, 7.1, 6.5, 7.4, 6.8, 7.0, 6.3, 7.2, 6.9, 6.6, 7.3, 6.1]

# Welch's t-test (default)
t, p = stats.ttest_ind(tumor, normal, equal_var=False)
print(f"Welch's t = {t:.4f}, p = {p:.2e}")

# Paired t-test
before = [245, 312, 198, 367, 289, 421, 156, 334, 278, 305]
after  = [180, 245, 165, 298, 220, 350, 132, 270, 210, 248]
t, p = stats.ttest_rel(before, after)
print(f"Paired t = {t:.4f}, p = {p:.6f}")

# Cohen's d (manual)
pooled_std = np.sqrt((np.std(tumor, ddof=1)**2 + np.std(normal, ddof=1)**2) / 2)
d = (np.mean(tumor) - np.mean(normal)) / pooled_std
print(f"Cohen's d = {d:.3f}")

R:

tumor  <- c(3.8, 4.5, 4.1, 3.2, 4.8, 3.9, 4.3, 5.1, 3.6, 4.0, 4.7, 3.5)
normal <- c(6.2, 7.1, 6.5, 7.4, 6.8, 7.0, 6.3, 7.2, 6.9, 6.6, 7.3, 6.1)

# Welch's t-test (default in R)
t.test(tumor, normal)

# Paired t-test
before <- c(245, 312, 198, 367, 289, 421, 156, 334, 278, 305)
after  <- c(180, 245, 165, 298, 220, 350, 132, 270, 210, 248)
t.test(before, after, paired = TRUE)

# Cohen's d
library(effsize)
cohen.d(tumor, normal)

Exercises

Exercise 1: Choose the Right t-Test

For each scenario, state which t-test variant is appropriate and why:

a) Comparing white blood cell counts between 20 patients with sepsis and 25 healthy volunteers b) Measuring gene expression in liver biopsies taken before and after drug treatment (same 15 patients) c) Testing whether mean read length from your sequencer matches the expected 150 bp

Exercise 2: Full t-Test Workflow

Hemoglobin levels (g/dL) in two groups:

• Anemia patients: [9.2, 8.8, 10.1, 9.5, 8.3, 9.7, 8.6, 9.0, 10.3, 8.9]
• Healthy controls: [13.5, 14.2, 12.8, 13.9, 14.5, 13.1, 14.0, 13.6, 12.9, 14.3]

let anemia  = [9.2, 8.8, 10.1, 9.5, 8.3, 9.7, 8.6, 9.0, 10.3, 8.9]
let healthy = [13.5, 14.2, 12.8, 13.9, 14.5, 13.1, 14.0, 13.6, 12.9, 14.3]

# TODO: 1. Check normality with qq_plot() on each group
# TODO: 2. Check equal variances by comparing variance() per group
# TODO: 3. Run the appropriate t-test with ttest()
# TODO: 4. Compute Cohen's d inline: (mean(a)-mean(b)) / sqrt((variance(a)+variance(b))/2)
# TODO: 5. Create a boxplot
# TODO: 6. Interpret results in a biological context

Exercise 3: Paired vs Independent

Run both a paired and independent t-test on the tumor volume data below. Compare the p-values and explain why they differ.

let before = [245, 312, 198, 367, 289, 421, 156, 334, 278, 305]
let after  = [180, 245, 165, 298, 220, 350, 132, 270, 210, 248]

# TODO: Run ttest_paired and ttest
# TODO: Which gives a smaller p-value? Why?
# TODO: What does the paired test "remove" that the independent test cannot?

Exercise 4: When Assumptions Fail

The following data are highly skewed (as often seen in cytokine measurements):

let treatment = [2.1, 1.8, 45.2, 3.5, 2.9, 1.2, 38.7, 4.1, 2.3, 1.5]
let control   = [0.8, 0.5, 0.9, 0.3, 1.1, 0.7, 0.4, 0.6, 1.0, 0.2]

# TODO: Test normality with qq_plot()
# TODO: Run the t-test with ttest() anyway — what does it say?
# TODO: Try log-transforming the data and re-testing
# TODO: Preview: tomorrow we'll learn non-parametric alternatives

Key Takeaways

• The t-test compares two group means, accounting for variability and sample size
• Welch’s t-test (unequal variances) should be the default — it is robust even when variances are equal
• Paired t-tests are more powerful when observations are naturally matched (same patient, same timepoint)
• Always check assumptions: Shapiro-Wilk for normality, Levene’s for equal variances, QQ plots for visual inspection
• Cohen’s d quantifies effect size independently of sample size: 0.2 = small, 0.5 = medium, 0.8 = large
• A significant t-test with a small Cohen’s d may not be biologically meaningful
• A non-significant t-test with a large Cohen’s d suggests you need more samples

What’s Next

What happens when your data violate the normality assumption? Cytokine levels, bacterial abundances, and many other biological measurements are wildly skewed. Tomorrow we introduce non-parametric tests — rank-based alternatives to the t-test that make no assumptions about the shape of your data distribution.

Day 9: When Normality Fails — Non-Parametric Tests

The Problem

Dr. Maria Gonzalez studies the gut microbiome in inflammatory bowel disease (IBD). She has 16S rRNA sequencing data from 15 IBD patients and 15 healthy controls, measuring the relative abundance of Faecalibacterium prausnitzii, a key anti-inflammatory bacterium. Looking at the data, she sees a mess: most values cluster near zero, a few patients have moderate levels, and one healthy individual has an enormous abundance of 45%. The histogram looks nothing like a bell curve — it is right-skewed with a long tail.

She runs a Shapiro-Wilk test on each group: both return p < 0.001, firmly rejecting normality. The t-test assumes normally distributed data. With data this skewed, the t-test’s p-value could be wildly inaccurate — too liberal or too conservative, depending on the specific pattern. She needs tests that work without any assumptions about the shape of the distribution.

These are non-parametric tests: methods that operate on the ranks of data rather than the raw values, making them robust to skewness, outliers, and any distributional shape.

What Are Non-Parametric Tests?

Imagine you are judging a cooking competition. A parametric judge scores each dish on a precise 1-100 scale and compares average scores. A non-parametric judge simply ranks the dishes from best to worst — first place, second place, third place. The ranking approach is less precise when scores are reliable, but it is far more robust when one judge has an eccentric scoring system.

Non-parametric tests replace raw data values with their ranks (1st smallest, 2nd smallest, …) and then analyze the ranks. This has powerful consequences:

Key insight: Non-parametric tests are not “worse” versions of parametric tests. They are the correct choice when distributional assumptions are violated. Using a t-test on heavily skewed data is like measuring temperature with a ruler — you might get a number, but it doesn’t mean anything.

When to Choose Non-Parametric

Use non-parametric tests when:

• Shapiro-Wilk rejects normality (p < 0.05) and sample size is small
• Data are ordinal (pain scale 1-10, tumor grade I-IV)
• Data have heavy outliers that cannot be removed
• Sample sizes are very small (n < 10 per group)
• Data are bounded or have floor/ceiling effects (many zeros)

The Rank Transformation

The foundation of all non-parametric tests is replacing values with ranks:

Notice: the outlier (45%) gets rank 6 — just one rank above 3.2%. Its extreme value no longer dominates the analysis.

Wilcoxon Rank-Sum Test (Mann-Whitney U)

The non-parametric counterpart of the independent two-sample t-test.

Procedure:

1. Combine all observations and rank them 1 through N
2. Sum the ranks in each group separately
3. If one group consistently has higher ranks, the rank sum will be extreme
4. Compare to the expected rank sum under H0 (no difference)

H0: The two groups have identical distributions H1: One group tends to have larger values

Common pitfall: The Wilcoxon rank-sum and Mann-Whitney U are the same test, just computed differently. U = W - n1(n1+1)/2. Different software uses different names, but the p-value is identical.

Wilcoxon Signed-Rank Test

The non-parametric counterpart of the paired t-test.

Procedure:

1. Compute the difference for each pair
2. Rank the absolute differences (ignoring zeros)
3. Sum ranks of positive differences (W+) and negative differences (W-)
4. If the treatment consistently increases (or decreases), one sum will dominate

Sign Test

Even simpler than Wilcoxon signed-rank — only considers the direction of differences, not their magnitude.

Procedure:

1. For each pair, note whether the difference is positive, negative, or zero
2. Count positives and negatives (discard zeros)
3. Under H0, positives and negatives should be equally likely (binomial test with p = 0.5)

The sign test has less power than Wilcoxon signed-rank but makes even fewer assumptions.

Kruskal-Wallis Test

The non-parametric counterpart of one-way ANOVA, for comparing three or more groups.

H0: All groups have the same distribution H1: At least one group differs

If significant, follow up with pairwise Wilcoxon tests (with multiple testing correction).

Kolmogorov-Smirnov (KS) Test

Compares two entire distributions, not just their centers. Detects differences in shape, spread, or location.

H0: The two samples come from the same distribution H1: The distributions differ in any way

Clinical relevance: The KS test is useful when you suspect groups differ not just in average abundance, but in the entire pattern of their distribution — for example, one group might be bimodal while the other is unimodal.

Decision Guide: Parametric vs Non-Parametric

Non-Parametric Tests in BioLang

Mann-Whitney U: Microbiome Abundance

# F. prausnitzii relative abundance (%) in IBD vs healthy
let ibd = [0.1, 0.3, 0.0, 0.8, 0.2, 0.0, 1.5, 0.4, 0.1, 0.0,
           3.2, 0.5, 0.1, 0.7, 0.0]
let healthy = [2.1, 5.4, 1.8, 8.2, 3.5, 12.1, 4.7, 6.3, 2.9, 45.0,
               7.1, 3.8, 9.5, 4.2, 6.8]

# First, demonstrate why t-test is inappropriate
# Check normality visually — both distributions are right-skewed
qq_plot(ibd, {title: "QQ Plot: IBD"})
qq_plot(healthy, {title: "QQ Plot: Healthy"})
print("Both groups are heavily skewed — normality violated!\n")

# Mann-Whitney U test (non-parametric)
let result = wilcoxon(ibd, healthy)
print("=== Mann-Whitney U Test ===")
print("U statistic: {result.statistic:.1}")
print("p-value: {result.p_value:.2e}")

# Compare to (inappropriate) t-test
let t_result = ttest(ibd, healthy)
print("\n(Inappropriate) Welch's t-test p-value: {t_result.p_value:.2e}")
print("Mann-Whitney p-value: {result.p_value:.2e}")
print("Results may differ substantially with skewed data")

# Visualize the skewed distributions
let bp_table = table({"IBD": ibd, "Healthy": healthy})
boxplot(bp_table, {title: "F. prausnitzii Abundance"})

Wilcoxon Signed-Rank: Paired Treatment Data

# Inflammatory cytokine IL-6 (pg/mL) before and after anti-TNF therapy
# Same 12 patients measured twice — highly skewed cytokine data
let before = [245, 18, 892, 45, 32, 1250, 67, 128, 15, 543, 78, 2100]
let after  = [120, 12, 340, 22, 28, 450,  35,  65, 10, 210, 42,  890]

# Normality check on differences
let diffs = zip(before, after) |> map(|p| p[0] - p[1])
qq_plot(diffs, {title: "QQ Plot: Paired Differences"})
print("Differences are non-normal -> use Wilcoxon signed-rank\n")

# Wilcoxon signed-rank test
let result = wilcoxon(before, after)
print("=== Wilcoxon Signed-Rank Test ===")
print("V statistic: {result.statistic:.1}")
print("p-value: {result.p_value:.6}")
print("All 12 patients showed reduction in IL-6")

# For comparison: the sign test via dbinom (even more robust, less powerful)
# Count how many differences are positive
let n_pos = diffs |> filter(|d| d > 0) |> len()
let n_nonzero = diffs |> filter(|d| d != 0) |> len()
# Under H0, n_pos ~ Binomial(n_nonzero, 0.5)
let sign_p = 0.0
for k in range(n_pos, n_nonzero + 1) {
    sign_p = sign_p + dbinom(k, n_nonzero, 0.5)
}
let sign_p = 2.0 * min(sign_p, 1.0 - sign_p)  # two-tailed
print("\nSign test p-value: {sign_p:.6}")

Kruskal-Wallis: Multiple Body Sites

# Bacterial diversity (Shannon index) across three gut regions
let ileum   = [1.2, 0.8, 1.5, 0.3, 2.1, 0.9, 1.4, 0.6, 1.8, 0.4]
let cecum   = [2.5, 3.1, 2.8, 2.2, 3.4, 2.7, 3.0, 2.3, 2.9, 3.2]
let rectum  = [3.8, 4.2, 3.5, 4.5, 3.9, 4.1, 3.6, 4.3, 3.7, 4.0]

# Kruskal-Wallis: use anova() on rank-transformed data
let result = anova([ileum, cecum, rectum])
print("=== Kruskal-Wallis Test: Diversity Across Body Sites ===")
print("H statistic: {result.statistic:.4}")
print("p-value: {result.p_value:.2e}")
print("Degrees of freedom: {result.df}")

if result.p_value < 0.05 {
  print("\nAt least one body site differs. Running pairwise comparisons...")

  let p1 = wilcoxon(ileum, cecum).p_value
  let p2 = wilcoxon(ileum, rectum).p_value
  let p3 = wilcoxon(cecum, rectum).p_value

  # Bonferroni correction for 3 comparisons
  let adjusted = p_adjust([p1, p2, p3], "bonferroni")
  print("Ileum vs Cecum:  p = {adjusted[0]:.4}")
  print("Ileum vs Rectum: p = {adjusted[1]:.4}")
  print("Cecum vs Rectum: p = {adjusted[2]:.4}")
}

let bp_table = table({"Ileum": ileum, "Cecum": cecum, "Rectum": rectum})
boxplot(bp_table, {title: "Microbial Diversity by Gut Region"})

KS Test: Comparing Distributions

# Do tumor suppressor genes and oncogenes have different
# expression distributions (not just different means)?
let tumor_suppressors = [2.1, 3.4, 1.8, 4.2, 2.9, 3.1, 2.5, 3.8, 1.5, 4.0,
                         2.7, 3.3, 2.0, 3.6, 2.3, 3.9, 1.9, 4.1, 2.6, 3.5]
let oncogenes = [5.2, 8.1, 6.3, 12.4, 7.5, 5.8, 9.2, 6.7, 11.3, 7.0,
                 5.5, 8.8, 6.1, 10.5, 7.3, 5.9, 9.7, 6.5, 11.8, 7.8]

let result = ks_test(tumor_suppressors, oncogenes)
print("=== KS Test: Expression Distributions ===")
print("D statistic: {result.statistic:.4}")
print("p-value: {result.p_value:.2e}")
print("Maximum distance between cumulative distributions: {result.statistic:.4}")

histogram([tumor_suppressors, oncogenes], {labels: ["Tumor Suppressors", "Oncogenes"], title: "Expression Distributions by Gene Class", x_label: "Expression (log2-CPM)", bins: 12})

Comparing t-Test vs Wilcoxon on the Same Data

set_seed(42)
# Demonstrate: with normal data, both tests agree
# With skewed data, they can disagree

print("=== Normal Data: Both Tests Agree ===")
let norm_a = rnorm(20, 5.0, 1.0)
let norm_b = rnorm(20, 6.0, 1.0)
let t_p = ttest(norm_a, norm_b).p_value
let w_p = wilcoxon(norm_a, norm_b).p_value
print("t-test p = {t_p:.4}, Mann-Whitney p = {w_p:.4}")

print("\n=== Skewed Data with Outlier: Tests May Disagree ===")
let skew_a = [1.2, 1.5, 1.8, 1.1, 1.4, 1.6, 1.3, 1.7, 1.9, 50.0]
let skew_b = [2.1, 2.3, 2.5, 2.0, 2.4, 2.2, 2.6, 2.1, 2.3, 2.5]
let t_p2 = ttest(skew_a, skew_b).p_value
let w_p2 = wilcoxon(skew_a, skew_b).p_value
print("t-test p = {t_p2:.4}, Mann-Whitney p = {w_p2:.4}")
print("The outlier inflates the t-test mean, masking the real pattern")

Python:

from scipy import stats

ibd = [0.1, 0.3, 0.0, 0.8, 0.2, 0.0, 1.5, 0.4, 0.1, 0.0,
       3.2, 0.5, 0.1, 0.7, 0.0]
healthy = [2.1, 5.4, 1.8, 8.2, 3.5, 12.1, 4.7, 6.3, 2.9, 45.0,
           7.1, 3.8, 9.5, 4.2, 6.8]

# Mann-Whitney U
u, p = stats.mannwhitneyu(ibd, healthy, alternative='two-sided')
print(f"U = {u}, p = {p:.2e}")

# Wilcoxon signed-rank (paired)
before = [245, 18, 892, 45, 32, 1250, 67, 128, 15, 543, 78, 2100]
after  = [120, 12, 340, 22, 28, 450,  35,  65, 10, 210, 42,  890]
w, p = stats.wilcoxon(before, after)
print(f"W = {w}, p = {p:.6f}")

# Kruskal-Wallis
ileum  = [1.2, 0.8, 1.5, 0.3, 2.1, 0.9, 1.4, 0.6, 1.8, 0.4]
cecum  = [2.5, 3.1, 2.8, 2.2, 3.4, 2.7, 3.0, 2.3, 2.9, 3.2]
rectum = [3.8, 4.2, 3.5, 4.5, 3.9, 4.1, 3.6, 4.3, 3.7, 4.0]
h, p = stats.kruskal(ileum, cecum, rectum)
print(f"H = {h:.4f}, p = {p:.2e}")

R:

ibd <- c(0.1, 0.3, 0.0, 0.8, 0.2, 0.0, 1.5, 0.4, 0.1, 0.0,
         3.2, 0.5, 0.1, 0.7, 0.0)
healthy <- c(2.1, 5.4, 1.8, 8.2, 3.5, 12.1, 4.7, 6.3, 2.9, 45.0,
             7.1, 3.8, 9.5, 4.2, 6.8)

wilcox.test(ibd, healthy)           # Mann-Whitney
wilcox.test(before, after, paired = TRUE)  # Wilcoxon signed-rank
kruskal.test(list(ileum, cecum, rectum))   # Kruskal-Wallis
ks.test(tumor_suppressors, oncogenes)      # KS test

Exercises

Exercise 1: Choose the Right Test

For each dataset, decide whether a parametric or non-parametric test is more appropriate:

a) Pain scores (0-10 scale) in drug vs placebo groups b) Blood pressure measurements in 30 patients (continuous, approximately normal) c) Number of bacterial colonies per plate (many zeros, some very high counts) d) Survival time in days (typically right-skewed)

Exercise 2: Microbiome Comparison

Two diets were compared for their effect on Bacteroides abundance:

let high_fiber = [8.2, 12.5, 6.3, 15.1, 9.8, 22.4, 7.5, 11.2, 14.8, 5.9]
let low_fiber  = [1.2, 0.5, 3.1, 0.8, 2.4, 0.3, 1.8, 0.9, 2.7, 0.6]

# TODO: Test normality of each group
# TODO: Run Mann-Whitney U test
# TODO: Also run a t-test and compare results
# TODO: Create a boxplot

Exercise 3: Multiple Body Sites with Post-Hoc

OTU richness from four body sites (oral, gut, skin, vaginal). Run Kruskal-Wallis and, if significant, perform all pairwise comparisons with Bonferroni correction.

let oral    = [120, 95, 145, 110, 88, 132, 105, 98, 140, 115]
let gut     = [350, 420, 280, 390, 310, 445, 360, 295, 410, 380]
let skin    = [180, 210, 165, 195, 220, 175, 200, 185, 230, 190]
let vaginal = [45, 30, 55, 38, 25, 50, 42, 35, 48, 28]

# TODO: Kruskal-Wallis test
# TODO: If significant, pairwise Mann-Whitney with Bonferroni correction
# TODO: Which sites differ from which?

Exercise 4: The Power Trade-Off

Generate 1000 simulations where both groups are truly normal with different means. Compare how often the t-test and Mann-Whitney detect the difference (power). Then repeat with skewed data (e.g., exponential).

# TODO: Simulate normal data, compare t-test vs Mann-Whitney power
# TODO: Simulate skewed data, compare again
# TODO: Which test wins in each scenario?

Key Takeaways

• Non-parametric tests use ranks instead of raw values, making them robust to skewness and outliers
• The Mann-Whitney U (Wilcoxon rank-sum) is the non-parametric alternative to the independent t-test
• The Wilcoxon signed-rank test is the non-parametric alternative to the paired t-test
• The Kruskal-Wallis test extends to three or more groups (non-parametric ANOVA)
• The KS test compares entire distributions, not just central tendency
• Non-parametric tests have about 95% of the power of parametric tests when data ARE normal, but are far more reliable when data are NOT normal
• Microbiome data, cytokine levels, survival times, and ordinal scales almost always require non-parametric methods
• Always check normality first (Shapiro-Wilk, QQ plots) — let the data guide your choice of test

What’s Next

So far we have compared two groups. But what if you have three, four, or ten groups — different drug doses, tissue types, or experimental conditions? Running all pairwise t-tests inflates false positives dramatically. Tomorrow we introduce ANOVA, the principled way to compare many groups simultaneously, along with post-hoc tests that identify which groups differ.

Day 10: Comparing Many Groups — ANOVA and Beyond

The Problem

Dr. James Park’s oncology team is testing a new targeted therapy at four dose levels: 0 mg (placebo), 25 mg, 50 mg, and 100 mg. Each group has 8 mice, and after 4 weeks they measure tumor volume in cubic millimeters. The team lead suggests: “Just do t-tests between all pairs of doses — that’s 6 comparisons, no big deal.”

But Dr. Park knows this is a trap. With 6 independent tests at alpha = 0.05, the probability of at least one false positive is not 5% — it is 1 - (0.95)^6 = 26.5%. Run 10 comparisons and it climbs to 40%. With 20,000 genes, the problem becomes catastrophic (we will tackle that on Day 12). The solution for comparing several groups at once is Analysis of Variance — ANOVA — which tests all groups simultaneously in a single, principled framework.

ANOVA has been the workhorse of experimental biology for nearly a century. Every drug dose-response study, every multi-tissue gene expression comparison, and every agricultural field trial relies on it. Today you will learn why it works, when it fails, and what to do after you get a significant result.

What Is ANOVA?

Imagine a classroom of students from four different schools. ANOVA asks: “Is the variation in test scores between schools larger than what we would expect given the variation within each school?”

If all four schools have similar average scores, the between-school variation will be small relative to within-school variation. If one school’s students consistently outscore the others, the between-school variation will dominate.

ANOVA decomposes total variation into two sources:

SS_total = SS_between + SS_within

The F-Statistic

The F-statistic is the ratio of between-group variance to within-group variance:

F = MS_between / MS_within

Where MS (mean square) = SS / df.

k = number of groups, N = total sample size.

• F near 1: Between-group variation is similar to within-group noise. No evidence of differences.
• F much greater than 1: Between-group variation exceeds what noise alone would produce. At least one group differs.

Key insight: ANOVA’s null hypothesis is that ALL group means are equal. A significant F-statistic tells you “at least one group differs” but does NOT tell you which one(s). You need post-hoc tests for that.

The Family-Wise Error Rate Problem

Why not just do multiple t-tests?

ANOVA controls this by testing all groups in a single hypothesis test.

Assumptions of One-Way ANOVA

1. Independence: Observations are independent within and between groups
2. Normality: Data within each group are approximately normally distributed
3. Homoscedasticity: All groups have equal variances (check with Levene’s test)

Common pitfall: ANOVA is robust to mild violations of normality when group sizes are equal and n > 10 per group. But it is sensitive to unequal variances, especially with unequal group sizes. When Levene’s test is significant, consider Welch’s ANOVA or the Kruskal-Wallis alternative.

Post-Hoc Tests: Which Groups Differ?

Tukey’s Honestly Significant Difference (HSD)

The gold standard post-hoc test. Compares all pairs of group means while controlling the family-wise error rate.

• Tests all k(k-1)/2 pairwise differences
• Provides adjusted p-values and confidence intervals
• Assumes equal variances and equal (or similar) group sizes

Other Post-Hoc Options

Effect Size: Eta-Squared

Just as Cohen’s d quantifies the effect for two groups, eta-squared quantifies it for ANOVA:

eta-squared = SS_between / SS_total

Non-Parametric Alternatives

ANOVA in BioLang

One-Way ANOVA: Dose-Response

# Tumor volume (mm^3) at 4 dose levels
let placebo = [485, 512, 468, 530, 495, 478, 521, 503]
let dose_25 = [420, 445, 398, 461, 432, 410, 452, 438]
let dose_50 = [310, 335, 288, 352, 321, 298, 345, 328]
let dose_100 = [180, 210, 165, 225, 195, 172, 218, 198]

# One-way ANOVA
let result = anova([placebo, dose_25, dose_50, dose_100])
print("=== One-Way ANOVA: Tumor Volume by Dose ===")
print("F-statistic: {result.statistic:.4}")
print("p-value: {result.p_value:.2e}")
print("df between: {result.df_between}, df within: {result.df_within}")

# Effect size: eta-squared = SS_between / SS_total (inline from anova output)
let eta2 = result.ss_between / (result.ss_between + result.ss_within)
print("Eta-squared: {eta2:.4} ({eta2*100:.1}% of variance explained)")

# Check assumptions: compare variances per group
print("\nVariances: placebo={variance(placebo):.1}, 25mg={variance(dose_25):.1}, 50mg={variance(dose_50):.1}, 100mg={variance(dose_100):.1}")

Tukey HSD Post-Hoc

let placebo  = [485, 512, 468, 530, 495, 478, 521, 503]
let dose_25  = [420, 445, 398, 461, 432, 410, 452, 438]
let dose_50  = [310, 335, 288, 352, 321, 298, 345, 328]
let dose_100 = [180, 210, 165, 225, 195, 172, 218, 198]

# Tukey HSD: pairwise ttest() + p_adjust()
let groups = [placebo, dose_25, dose_50, dose_100]
let labels = ["Placebo", "25mg", "50mg", "100mg"]
let pairwise_p = []
let pair_labels = []
for i in 0..len(groups) {
  for j in (i+1)..len(groups) {
    let r = ttest(groups[i], groups[j])
    pairwise_p = append(pairwise_p, r.p_value)
    pair_labels = append(pair_labels, "{labels[i]} vs {labels[j]}")
  }
}
let adj_p = p_adjust(pairwise_p, "bonferroni")

print("=== Pairwise t-tests (Bonferroni-adjusted) ===")
print("Comparison          | Diff     | p-adj")
print("--------------------|----------|--------")
for k in 0..len(pair_labels) {
  print("{pair_labels[k]:<20}| {adj_p[k]:.4}")
}

Grouped Boxplot Visualization

let groups = {
  "Placebo": [485, 512, 468, 530, 495, 478, 521, 503],
  "25 mg":   [420, 445, 398, 461, 432, 410, 452, 438],
  "50 mg":   [310, 335, 288, 352, 321, 298, 345, 328],
  "100 mg":  [180, 210, 165, 225, 195, 172, 218, 198]
}

boxplot(groups, {title: "Tumor Volume by Treatment Dose", y_label: "Tumor Volume (mm^3)", x_label: "Dose Group", show_points: true})

Kruskal-Wallis: When ANOVA Assumptions Fail

# Cytokine levels across three disease stages — heavily skewed
let stage_I   = [2.1, 1.5, 3.8, 0.9, 12.5, 1.8, 2.4, 0.7, 4.2, 1.1]
let stage_II  = [8.5, 15.2, 5.3, 22.1, 9.8, 45.0, 7.2, 12.8, 6.1, 18.5]
let stage_III = [35.2, 88.1, 42.5, 120.0, 55.3, 78.9, 95.2, 48.7, 110.5, 65.8]

# Check normality
# Visual normality check — all stages are right-skewed
qq_plot(stage_I, {title: "QQ: Stage I"})
qq_plot(stage_II, {title: "QQ: Stage II"})
qq_plot(stage_III, {title: "QQ: Stage III"})
print("Normality violated -> use Kruskal-Wallis (anova on ranks)\n")

let result = anova([stage_I, stage_II, stage_III])
print("=== Kruskal-Wallis Test ===")
print("H statistic: {result.statistic:.4}")
print("p-value: {result.p_value:.2e}")

# Pairwise follow-up with Bonferroni correction
if result.p_value < 0.05 {
  let p12 = wilcoxon(stage_I, stage_II).p_value
  let p13 = wilcoxon(stage_I, stage_III).p_value
  let p23 = wilcoxon(stage_II, stage_III).p_value
  let adj = p_adjust([p12, p13, p23], "bonferroni")

  print("\nPairwise Mann-Whitney (Bonferroni-adjusted):")
  print("  Stage I vs II:   p = {adj[0]:.4}")
  print("  Stage I vs III:  p = {adj[1]:.4}")
  print("  Stage II vs III: p = {adj[2]:.4}")
}

let bp_table = table({"Stage I": stage_I, "Stage II": stage_II, "Stage III": stage_III})
boxplot(bp_table, {title: "IL-6 Levels by Disease Stage"})

Friedman Test: Repeated Measures

# Pain scores (1-10) for 8 patients under 3 analgesics
# Same patients tested with each drug (crossover design)
let drug_a = [7, 5, 8, 6, 4, 7, 5, 6]
let drug_b = [4, 3, 5, 4, 2, 5, 3, 4]
let drug_c = [3, 2, 4, 3, 1, 3, 2, 3]

# Friedman test: anova() on ranks for repeated measures
let result = anova([drug_a, drug_b, drug_c])
print("=== Friedman Test: Pain Scores Across Analgesics ===")
print("Chi-squared: {result.statistic:.4}")
print("p-value: {result.p_value:.6}")

if result.p_value < 0.05 {
  print("At least one analgesic differs in pain reduction")
  print("Medians: Drug A={median(drug_a)}, Drug B={median(drug_b)}, Drug C={median(drug_c)}")
}

Complete Workflow: Gene Expression Across Tissues

# FOXP3 expression across immune cell types
let t_reg  = [8.5, 9.2, 8.8, 9.5, 8.1, 9.0, 8.7, 9.3]
let t_eff  = [3.2, 3.8, 3.5, 4.1, 2.9, 3.6, 3.3, 3.9]
let b_cell = [1.5, 1.8, 1.2, 2.0, 1.4, 1.7, 1.3, 1.9]
let nk     = [0.8, 1.1, 0.6, 1.3, 0.9, 1.0, 0.7, 1.2]

# Step 1: Check assumptions
print("=== Assumption Checks ===")
# Compare variances across groups
print("Variances: T-reg={variance(t_reg):.3}, T-eff={variance(t_eff):.3}, B cell={variance(b_cell):.3}, NK={variance(nk):.3}")
# Visual normality check
for name, data in [["T-reg", t_reg], ["T-eff", t_eff], ["B cell", b_cell], ["NK", nk]] {
  qq_plot(data, {title: "QQ Plot: {name}"})
}

# Step 2: ANOVA
let result = anova([t_reg, t_eff, b_cell, nk])
print("\n=== One-Way ANOVA ===")
print("F = {result.statistic:.2}, p = {result.p_value:.2e}")

# Step 3: Effect size (eta-squared from anova output)
let eta2 = result.ss_between / (result.ss_between + result.ss_within)
print("Eta-squared = {eta2:.3} (cell type explains {eta2*100:.1}% of variance)")

# Step 4: Post-hoc — pairwise ttest() + p_adjust()
let cell_groups = [t_reg, t_eff, b_cell, nk]
let cell_labels = ["T-reg", "T-eff", "B cell", "NK"]
let pw_pvals = []
let pw_labels = []
for i in 0..len(cell_groups) {
  for j in (i+1)..len(cell_groups) {
    pw_pvals = append(pw_pvals, ttest(cell_groups[i], cell_groups[j]).p_value)
    pw_labels = append(pw_labels, "{cell_labels[i]} vs {cell_labels[j]}")
  }
}
let pw_adj = p_adjust(pw_pvals, "bonferroni")

print("\n=== Pairwise t-tests (Bonferroni) ===")
for k in 0..len(pw_labels) {
  let sig = if pw_adj[k] < 0.001 then "***" else if pw_adj[k] < 0.01 then "**" else if pw_adj[k] < 0.05 then "*" else "ns"
  print("{pw_labels[k]:<20} p={pw_adj[k]:.4} {sig}")
}

# Step 5: Visualize
let bp_table = table({"T-reg": t_reg, "T-eff": t_eff, "B cell": b_cell, "NK": nk})
boxplot(bp_table, {title: "FOXP3 Expression Across Immune Cell Types", show_points: true})

Python:

from scipy import stats
import scikit_posthocs as sp

placebo  = [485, 512, 468, 530, 495, 478, 521, 503]
dose_25  = [420, 445, 398, 461, 432, 410, 452, 438]
dose_50  = [310, 335, 288, 352, 321, 298, 345, 328]
dose_100 = [180, 210, 165, 225, 195, 172, 218, 198]

# One-way ANOVA
f, p = stats.f_oneway(placebo, dose_25, dose_50, dose_100)
print(f"F = {f:.4f}, p = {p:.2e}")

# Tukey HSD
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import numpy as np
data = placebo + dose_25 + dose_50 + dose_100
groups = ['P']*8 + ['25']*8 + ['50']*8 + ['100']*8
print(pairwise_tukeyhsd(data, groups))

# Kruskal-Wallis
h, p = stats.kruskal(placebo, dose_25, dose_50, dose_100)

# Friedman
stats.friedmanchisquare(drug_a, drug_b, drug_c)

R:

# One-way ANOVA
data <- data.frame(
  volume = c(485,512,468,530,495,478,521,503,
             420,445,398,461,432,410,452,438,
             310,335,288,352,321,298,345,328,
             180,210,165,225,195,172,218,198),
  dose = factor(rep(c("Placebo","25mg","50mg","100mg"), each=8))
)
result <- aov(volume ~ dose, data = data)
summary(result)
TukeyHSD(result)

# Eta-squared
library(effectsize)
eta_squared(result)

# Kruskal-Wallis
kruskal.test(volume ~ dose, data = data)

# Friedman
friedman.test(matrix(c(drug_a, drug_b, drug_c), ncol=3))

Exercises

Exercise 1: FWER Calculation

You have 5 experimental groups and want to compare all pairs. Calculate: (a) how many pairwise comparisons there are, (b) the probability of at least one false positive at alpha = 0.05, (c) what the Bonferroni-adjusted alpha would be.

Exercise 2: Full ANOVA Workflow

Three fertilizer treatments were tested on plant growth (cm):

let control   = [12.5, 14.2, 11.8, 13.5, 12.9, 14.8, 11.2, 13.1]
let fert_a    = [18.3, 20.1, 17.5, 19.8, 18.9, 21.2, 17.0, 19.5]
let fert_b    = [15.8, 17.2, 14.5, 16.9, 15.3, 17.8, 14.1, 16.5]

# TODO: Check normality and equal variances
# TODO: Run one-way ANOVA
# TODO: Compute eta-squared
# TODO: If significant, run Tukey HSD
# TODO: Create boxplot
# TODO: Interpret: which fertilizer is best?

Exercise 3: Parametric vs Non-Parametric ANOVA

Run both ANOVA and Kruskal-Wallis on the cytokine data. Compare p-values. Which is more appropriate?

let mild     = [5.2, 3.8, 8.1, 2.5, 12.0, 4.3, 6.7, 1.9]
let moderate = [25.1, 45.0, 18.3, 52.8, 30.5, 15.2, 38.7, 22.4]
let severe   = [120, 250, 85, 310, 180, 95, 275, 145]

# TODO: Check normality
# TODO: Run both ANOVA and Kruskal-Wallis
# TODO: Which test is more appropriate for this data? Why?

Exercise 4: Repeated Measures Design

Five patients had their blood pressure measured under three conditions (rest, mild exercise, intense exercise):

let rest     = [120, 135, 118, 142, 125]
let mild_ex  = [130, 148, 128, 155, 138]
let intense  = [155, 172, 148, 180, 162]

# TODO: Use Friedman test (non-parametric repeated measures)
# TODO: If significant, perform pairwise Wilcoxon signed-rank tests
# TODO: Apply Bonferroni correction to the pairwise p-values

Key Takeaways

• Multiple t-tests inflate the false positive rate — the family-wise error rate grows rapidly with the number of comparisons
• ANOVA tests whether any group differs from the others in a single F-test, controlling the overall error rate
• The F-statistic compares between-group variance to within-group variance: F much greater than 1 suggests real differences
• A significant ANOVA tells you “at least one group differs” — use Tukey HSD post-hoc to find which pairs differ
• Eta-squared measures effect size: the proportion of total variance explained by group membership
• Kruskal-Wallis is the non-parametric alternative when normality is violated
• Friedman test handles repeated measures designs non-parametrically
• Always check assumptions (normality, equal variances) before interpreting ANOVA results

What’s Next

Tomorrow we shift from continuous to categorical outcomes. When your data consist of counts in categories — genotypes, disease status, response/non-response — you need the chi-square test and Fisher’s exact test. These tools are essential for testing genetic associations, evaluating Hardy-Weinberg equilibrium, and computing odds ratios in case-control studies.

Day 11: Categorical Data — Chi-Square and Fisher’s Exact

The Problem

Dr. Elena Vasquez is an epidemiologist studying the genetics of Alzheimer’s disease. She genotypes a SNP near the APOE gene in 1,000 participants: 500 with Alzheimer’s and 500 age-matched controls. Among Alzheimer’s patients, 180 carry at least one copy of the risk allele. Among controls, 120 carry it. The numbers look different — 36% versus 24% — but these are proportions, not measurements on a continuous scale. She cannot compute a mean or standard deviation. She cannot run a t-test.

When your data are counts in categories — disease yes/no, genotype AA/AG/GG, response/non-response — you need tests designed for categorical data. The chi-square test and Fisher’s exact test are the workhorses. They also underpin some of the most important calculations in genetics: Hardy-Weinberg equilibrium, odds ratios for case-control studies, and allelic association tests.

This chapter covers the full toolkit for analyzing categorical data, from contingency tables to effect size measures like odds ratios and Cramer’s V.

What Are Categorical Data?

Categorical variables place observations into discrete groups rather than measuring them on a continuous scale.

The fundamental data structure for categorical analysis is the contingency table (also called a cross-tabulation):

Chi-Square Test of Independence

The chi-square test asks: “Are these two categorical variables independent, or is there an association?”

How It Works

1. Compute expected counts under independence: E = (row total x column total) / grand total
2. For each cell, compute (Observed - Expected)^2 / Expected
3. Sum across all cells to get the chi-square statistic
4. Compare to the chi-square distribution with df = (rows - 1)(cols - 1)

For the Alzheimer’s example:

The chi-square statistic measures how far the observed counts deviate from what independence predicts.

Assumptions

• Observations are independent
• Expected count in each cell is at least 5 (if not, use Fisher’s exact test)
• Sample is reasonably large

Common pitfall: The chi-square test requires expected counts of at least 5 in each cell, not observed counts. Check expected counts before interpreting results. When they are too small, Fisher’s exact test is the safe alternative.

Chi-Square Goodness of Fit: Hardy-Weinberg

The chi-square test can also compare observed frequencies to a theoretical distribution. A critical application in genetics is testing Hardy-Weinberg Equilibrium (HWE).

For a biallelic locus with allele frequencies p (major) and q (minor), HWE predicts:

• AA: p^2
• Ag: 2pq
• gg: q^2

If observed genotype counts deviate significantly from HWE expectations, it may indicate selection, population structure, genotyping error, or non-random mating.

Fisher’s Exact Test

When sample sizes are small or expected counts fall below 5, the chi-square approximation is unreliable. Fisher’s exact test computes the exact probability using the hypergeometric distribution.

Fisher’s exact test is computationally expensive for large tables but is the gold standard for 2x2 tables with small counts — common in rare variant studies and pilot experiments.

Clinical relevance: Regulatory agencies (FDA, EMA) often prefer Fisher’s exact test for safety analyses where adverse events are rare and sample sizes are modest.

McNemar’s Test: Paired Categorical Data

When observations are paired — the same patients tested before and after treatment, or the same samples tested with two diagnostic methods — McNemar’s test is the correct choice.

McNemar’s test focuses on the discordant pairs (b and c):

chi-square = (b - c)^2 / (b + c)

Measures of Association

Odds Ratio (OR)

The odds ratio quantifies the strength of association in a 2x2 table:

OR = (a x d) / (b x c)

Relative Risk (RR)

In prospective studies (cohorts), relative risk is often preferred:

RR = [a / (a+b)] / [c / (c+d)]

Key insight: Odds ratios approximate relative risk only when the outcome is rare (< 10%). For common outcomes, OR overestimates RR. Case-control studies can only estimate OR, not RR.

Cramer’s V

A measure of association strength for tables larger than 2x2:

V = sqrt(chi-square / (n x min(r-1, c-1)))

V ranges from 0 (no association) to 1 (perfect association).

Proportion Test

Tests whether an observed proportion differs from an expected value, or whether two proportions differ from each other. Useful for comparing mutation frequencies, response rates, or allele frequencies between populations.

Categorical Tests in BioLang

Chi-Square Test: SNP-Disease Association

# Observed genotype counts near APOE
# Rows: Alzheimer's, Control
# Columns: Risk allele present, Risk allele absent
let observed = [[180, 320], [120, 380]]

let result = chi_square(observed)
print("=== Chi-Square Test of Independence ===")
print("Chi-square statistic: {result.statistic:.4}")
print("p-value: {result.p_value:.2e}")
print("Degrees of freedom: {result.df}")

# Effect size: Cramer's V from chi-square output
let n_total = 180 + 320 + 120 + 380
let v = sqrt(result.statistic / (n_total * min(2 - 1, 2 - 1)))
print("Cramer's V: {v:.4}")

# Odds ratio (inline: (a*d) / (b*c))
let a = 180
let b = 320
let c = 120
let d = 380
let or_val = (a * d) / (b * c)
print("\nOdds Ratio: {or_val:.3}")

# Relative risk (inline)
let rr = (a / (a + b)) / (c / (c + d))
print("Relative Risk: {rr:.3}")

if result.p_value < 0.05 {
  print("\nSignificant association between risk allele and Alzheimer's")
}

Fisher’s Exact Test: Rare Variant Study

# Rare loss-of-function variant in 200 cases and 200 controls
# Very small expected counts -> Fisher's exact test
let observed = [[8, 192], [2, 198]]

print("=== Fisher's Exact Test: Rare Variant ===")
print("Observed: 8/200 cases vs 2/200 controls carry the variant\n")

# Chi-square would be unreliable here
let chi_result = chi_square(observed)
print("Chi-square p-value: {chi_result.p_value:.4} (unreliable — low expected counts)")

# Fisher's exact is the correct choice
let fisher_result = fisher_exact(observed)
print("Fisher's exact p-value: {fisher_result.p_value:.4}")

# Odds ratio (inline)
let or_val = (8 * 198) / (192 * 2)
print("Odds Ratio: {or_val:.2}")

# Note: CI includes 1.0, so despite the apparent 4x difference,
# the sample size is too small for significance

Hardy-Weinberg Equilibrium Test

# Genotype counts at a SNP locus
# Observed: AA=280, AG=430, GG=290 (total = 1000)
let obs_AA = 280
let obs_AG = 430
let obs_GG = 290
let n = obs_AA + obs_AG + obs_GG

# Estimate allele frequencies
let p = (2 * obs_AA + obs_AG) / (2 * n)  # freq of A
let q = 1.0 - p                           # freq of G

print("Allele frequencies: p(A) = {p:.4}, q(G) = {q:.4}")

# Expected counts under HWE
let exp_AA = p * p * n
let exp_AG = 2 * p * q * n
let exp_GG = q * q * n

print("\nGenotype     | Observed | Expected (HWE)")
print("-------------|----------|----------------")
print("AA           | {obs_AA:>8} | {exp_AA:>14.1}")
print("AG           | {obs_AG:>8} | {exp_AG:>14.1}")
print("GG           | {obs_GG:>8} | {exp_GG:>14.1}")

# Chi-square goodness of fit (df=1 for HWE with 2 alleles)
let chi_sq = (obs_AA - exp_AA)^2 / exp_AA +
             (obs_AG - exp_AG)^2 / exp_AG +
             (obs_GG - exp_GG)^2 / exp_GG

# Use chi_square with expected frequencies
let observed = [obs_AA, obs_AG, obs_GG]
let expected = [exp_AA, exp_AG, exp_GG]
let result = chi_square(observed, expected)

print("\nChi-square = {result.statistic:.4}")
print("p-value = {result.p_value:.4}")

if result.p_value > 0.05 {
  print("Genotype frequencies are consistent with Hardy-Weinberg Equilibrium")
} else {
  print("Significant deviation from HWE — investigate possible causes")
}

McNemar’s Test: Diagnostic Agreement

# Two diagnostic tests for TB applied to same 200 patients
# Test A (culture) vs Test B (PCR)
let table = [[85, 15], [10, 90]]
# 85: both positive, 90: both negative
# 15: A+/B-, 10: A-/B+

# McNemar's test: use chi_square() on discordant cells
let b_disc = 15  # A+/B-
let c_disc = 10  # A-/B+
let mcnemar_chi2 = (b_disc - c_disc) * (b_disc - c_disc) / (b_disc + c_disc)
let mcnemar_p = 1.0 - pnorm(sqrt(mcnemar_chi2), 0, 1) * 2.0  # approximate
# Or use chi_square on discordant cells
let result = chi_square([b_disc, c_disc], [12.5, 12.5])
print("=== McNemar's Test: Culture vs PCR ===")
print("Discordant pairs: A+/B- = 15, A-/B+ = 10")
print("Chi-square: {result.statistic:.4}")
print("p-value: {result.p_value:.4}")

if result.p_value > 0.05 {
  print("No significant difference between the two diagnostic tests")
} else {
  print("The tests have significantly different detection rates")
}

Proportion Test: Comparing Mutation Frequencies

# EGFR mutation frequency in two populations
# Asian cohort: 120/300 (40%), European cohort: 45/300 (15%)
# Two-proportion test via chi_square
let observed = [[120, 180], [45, 255]]
let result = chi_square(observed)

print("=== Two-Proportion Test: EGFR Mutation Frequency ===")
print("Asian: 120/300 = 40%")
print("European: 45/300 = 15%")
print("Chi-square: {result.statistic:.4}")
print("p-value: {result.p_value:.2e}")
print("Difference: 25 percentage points")

# Visualize
bar_chart(["Asian", "European"], [40.0, 15.0], {title: "EGFR Mutation Frequency by Population", y_label: "Mutation Frequency (%)"})

Complete Workflow: Multi-Allelic Association

# Three genotypes at a pharmacogenomics locus
# and three drug response categories
let observed = [
  [45, 30, 25],   # Poor metabolizer
  [35, 55, 60],   # Intermediate
  [20, 15, 15]    # Ultra-rapid
]
let row_labels = ["Poor", "Intermediate", "Ultra-rapid"]
let col_labels = ["AA", "AG", "GG"]

let result = chi_square(observed)
print("=== Chi-Square: Genotype vs Drug Response ===")
print("Chi-square: {result.statistic:.4}")
print("p-value: {result.p_value:.4}")
print("df: {result.df}")

# Cramer's V from chi-square output
let n_total = 45 + 30 + 25 + 35 + 55 + 60 + 20 + 15 + 15
let v = sqrt(result.statistic / (n_total * min(3 - 1, 3 - 1)))
print("Cramer's V: {v:.4} (effect size)")

# Display the contingency table
print("\n           | AA   | AG   | GG   | Total")
print("-----------|------|------|------|------")
for i in 0..3 {
  let total = observed[i][0] + observed[i][1] + observed[i][2]
  print("{row_labels[i]:<11}| {observed[i][0]:>4} | {observed[i][1]:>4} | {observed[i][2]:>4} | {total:>4}")
}

Python:

from scipy import stats
import numpy as np

# Chi-square test of independence
observed = np.array([[180, 320], [120, 380]])
chi2, p, dof, expected = stats.chi2_contingency(observed)
print(f"Chi-square = {chi2:.4f}, p = {p:.2e}")

# Fisher's exact test
odds, p = stats.fisher_exact([[8, 192], [2, 198]])
print(f"OR = {odds:.2f}, p = {p:.4f}")

# McNemar's test
from statsmodels.stats.contingency_tables import mcnemar
result = mcnemar(np.array([[85, 15], [10, 90]]), exact=False)
print(f"McNemar chi2 = {result.statistic:.4f}, p = {result.pvalue:.4f}")

# Proportion test
from statsmodels.stats.proportion import proportions_ztest
z, p = proportions_ztest([120, 45], [300, 300])
print(f"z = {z:.4f}, p = {p:.2e}")

R:

# Chi-square test
observed <- matrix(c(180, 120, 320, 380), nrow = 2)
chisq.test(observed)

# Fisher's exact test
fisher.test(matrix(c(8, 2, 192, 198), nrow = 2))

# McNemar's test
mcnemar.test(matrix(c(85, 10, 15, 90), nrow = 2))

# Odds ratio
library(epitools)
oddsratio(observed)

# Proportion test
prop.test(c(120, 45), c(300, 300))

Exercises

Exercise 1: SNP Association Study

A GWAS hit is validated in a replication cohort. Genotype counts:

let observed = [[150, 200, 50], [180, 160, 60]]

# TODO: Run chi-square test
# TODO: Compute Cramer's V
# TODO: Test HWE separately in cases and controls
# TODO: Interpret the results

Exercise 2: Fisher’s Exact on Rare Mutations

In a rare disease study: 5/50 patients carry a mutation vs 1/100 controls.

let table = [[5, 45], [1, 99]]

# TODO: Why is Fisher's exact preferred here?
# TODO: Compute p-value and odds ratio
# TODO: What does the wide CI on the OR tell you?

Exercise 3: McNemar’s for Treatment Response

A tumor is biopsied before and after chemotherapy. Response to an immunostaining marker:

let table = [[40, 25], [5, 30]]

# TODO: Run McNemar's test
# TODO: What do the discordant pairs tell you?
# TODO: Is the marker significantly altered by chemotherapy?

Exercise 4: Multi-Population Allele Comparison

Compare allele frequencies of a pharmacogenomics variant across three populations. Use chi-square and Cramer’s V, then create a bar chart of frequencies.

# Observed carrier counts out of 500 per population
let observed = [[210, 290], [150, 350], [85, 415]]
let populations = ["East Asian", "European", "African"]

# TODO: Chi-square test of independence
# TODO: Cramer's V
# TODO: Pairwise proportion tests with Bonferroni correction
# TODO: Bar chart of carrier frequencies

Key Takeaways

• The chi-square test evaluates whether two categorical variables are independent by comparing observed to expected counts
• Fisher’s exact test is preferred when expected cell counts are below 5, common with rare variants
• McNemar’s test handles paired categorical data (same subjects, two conditions)
• The odds ratio quantifies association strength; OR = 1 means no association
• Relative risk is preferred in cohort studies but cannot be computed from case-control designs
• Cramer’s V measures association strength for tables larger than 2x2
• Hardy-Weinberg equilibrium testing uses chi-square goodness of fit to check for genotyping artifacts or population structure
• The proportion test compares frequencies between populations or against expected values
• Always check expected cell counts before using chi-square — use Fisher’s exact when they are small

What’s Next

You now have a powerful toolkit for individual tests. But in genomics, we never run just one test — we run thousands or millions simultaneously. Testing 20,000 genes means 1,000 false positives at alpha = 0.05. Tomorrow we confront the multiple testing crisis head-on and learn the corrections that make genome-scale analysis possible, culminating in the volcano plot that has become the icon of differential expression analysis.

Day 12: The Multiple Testing Crisis — FDR and Correction

The Problem

Dr. Rachel Kim is a genomicist analyzing differential gene expression between tumor and normal tissue. She runs a Welch’s t-test on each of 20,000 genes and finds 1,200 with p < 0.05. Exciting — until she does the arithmetic: 20,000 genes multiplied by a 5% false positive rate equals 1,000 genes that would appear “significant” purely by chance, even if not a single gene were truly differentially expressed.

Her 1,200 “hits” likely contain about 1,000 false positives and perhaps 200 true findings. That is a false discovery rate of over 80%. If she publishes this list and a validation lab tries to confirm the top 50 genes, roughly 40 will fail to replicate. Her scientific reputation — and the lab’s funding — depends on solving this problem.

This is the multiple testing crisis, and it is the single most important statistical concept in genomics. Every differential expression study, every GWAS, every proteomics screen must address it. The solutions — Bonferroni correction, Benjamini-Hochberg FDR, and their relatives — are what make genome-scale analysis possible.

Making It Visceral: A Simulation

Before discussing theory, let us see the problem with our own eyes.

set_seed(42)
# Simulate 20,000 genes where NONE are truly different (complete null)

let p_values = []
for i in 1..20000 {
  # Both groups drawn from the same distribution — no real differences
  let group1 = rnorm(10, 0, 1)
  let group2 = rnorm(10, 0, 1)
  let result = ttest(group1, group2)
  p_values = append(p_values, result.p_value)
}

# Count "discoveries" at various thresholds
let sig_05 = p_values |> filter(|p| p < 0.05) |> len()
let sig_01 = p_values |> filter(|p| p < 0.01) |> len()
let sig_001 = p_values |> filter(|p| p < 0.001) |> len()

print("=== 20,000 Null Genes (No True Differences) ===")
print("'Significant' at p < 0.05:  {sig_05} (expected: ~1000)")
print("'Significant' at p < 0.01:  {sig_01} (expected: ~200)")
print("'Significant' at p < 0.001: {sig_001} (expected: ~20)")
print("\nEvery single one is a false positive!")

histogram(p_values, {title: "p-Value Distribution Under Complete Null", x_label: "p-value", bins: 50})

Under the null hypothesis, p-values are uniformly distributed between 0 and 1. Exactly 5% will fall below 0.05 by definition. With 20,000 tests, that is 1,000 false alarms.

Family-Wise Error Rate (FWER)

The family-wise error rate is the probability of making at least one Type I error across all tests:

FWER = 1 - (1 - alpha)^m

Where m is the number of tests.

Bonferroni Correction

The simplest and most conservative approach: divide alpha by the number of tests.

Adjusted alpha = alpha / m

For 20,000 tests at alpha = 0.05: adjusted alpha = 0.05 / 20,000 = 2.5 x 10^-6.

Equivalently, multiply each p-value by m and compare to the original alpha:

p_adjusted = min(p x m, 1.0)

Strengths and Weaknesses

Common pitfall: Bonferroni is often TOO conservative for genomics. With 20,000 correlated gene expression measurements, it throws out far too many true positives. It is appropriate when you need to be extremely cautious (drug safety) or when you have few tests.

Holm’s Step-Down Correction

Holm’s method is uniformly more powerful than Bonferroni while still controlling FWER.

Procedure:

1. Sort p-values from smallest to largest: p(1) <= p(2) <= … <= p(m)
2. For the i-th smallest p-value, compute adjusted p: p(i) x (m - i + 1)
3. Enforce monotonicity: each adjusted p must be >= the previous
4. Reject all hypotheses whose adjusted p < alpha

Holm’s method is always at least as powerful as Bonferroni, and sometimes substantially more so.

False Discovery Rate (FDR): The Breakthrough

In 1995, Benjamini and Hochberg introduced a paradigm shift. Instead of controlling the probability of any false positive (FWER), they controlled the expected proportion of false positives among rejected hypotheses.

FDR = E[V / R]

Where V = number of false positives and R = total rejections.

If you set FDR = 0.05, you accept that about 5% of your “discoveries” may be false — a much more practical threshold for genomics than guaranteeing zero false positives.

Benjamini-Hochberg (BH) Procedure

Procedure:

1. Sort p-values from smallest to largest: p(1) <= p(2) <= … <= p(m)
2. For each p(i), compute the BH threshold: (i / m) x q, where q is the desired FDR level
3. Find the largest i where p(i) <= (i / m) x q
4. Reject all hypotheses 1, 2, …, i

Equivalently, the BH-adjusted p-value (q-value) is:

q(i) = min(p(i) x m / i, 1.0) (with monotonicity enforced)

Why BH Changed Genomics

Key insight: BH-FDR is the standard for differential expression, GWAS, proteomics, and almost all high-throughput biology. When a paper reports “genes with FDR < 0.05,” they almost always mean Benjamini-Hochberg adjusted p-values.

Other Correction Methods

Hochberg’s Step-Up

Similar to Holm but slightly more powerful; assumes independence or positive dependence among tests.

Benjamini-Yekutieli (BY)

A conservative FDR procedure that works under any dependency structure between tests. Use when you have strong correlations (e.g., genes in the same pathway).

Choosing a Method

The Volcano Plot

The volcano plot is the most iconic visualization in differential expression analysis. It plots:

• x-axis: log2 fold change (effect size)
• y-axis: -log10(adjusted p-value) (statistical significance)

Genes in the upper corners are both statistically significant AND biologically meaningful — the ones you actually care about.

Multiple Testing Correction in BioLang

Simulating the Crisis and Applying Corrections

set_seed(42)
# Simulate 20,000 genes: 18,000 null + 2,000 truly differential

let p_values = []
let is_true = []  # Ground truth: 1 = truly differential, 0 = null
let fold_changes = []

for i in 1..20000 {
  let group1 = rnorm(10, 0, 1)
  if i <= 2000 {
    # True differential: shifted mean
    let shift = rnorm(1, 2.0, 0.5)[0]
    let group2 = rnorm(10, shift, 1)
    is_true = append(is_true, 1)
  } else {
    # Null gene: no difference
    let group2 = rnorm(10, 0, 1)
    is_true = append(is_true, 0)
  }
  let result = ttest(group1, group2)
  p_values = append(p_values, result.p_value)
  fold_changes = append(fold_changes, mean(group2) - mean(group1))
}

print("Total genes: {len(p_values)}")
print("True DE genes: {is_true |> filter(|x| x == 1) |> len()}")
print("Null genes: {is_true |> filter(|x| x == 0) |> len()}")

Applying All Correction Methods

# Apply multiple correction methods
let p_bonf = p_adjust(p_values, "bonferroni")
let p_holm = p_adjust(p_values, "holm")
let p_bh   = p_adjust(p_values, "BH")
let p_hoch = p_adjust(p_values, "hochberg")
let p_by   = p_adjust(p_values, "BY")

# Count discoveries at adjusted p < 0.05
let count_sig = |adj_p| adj_p |> filter(|p| p < 0.05) |> len()

print("=== Discoveries at Adjusted p < 0.05 ===")
print("Method       | Total Discoveries | True Pos | False Pos | FDR")
print("-------------|-------------------|----------|-----------|--------")

for method_name, adj_p in [
  ["Unadjusted ", p_values],
  ["Bonferroni ", p_bonf],
  ["Holm       ", p_holm],
  ["BH (FDR)   ", p_bh],
  ["Hochberg   ", p_hoch],
  ["BY         ", p_by]
] {
  let discoveries = []
  let true_pos = 0
  let false_pos = 0
  for j in 0..len(adj_p) {
    if adj_p[j] < 0.05 {
      discoveries = append(discoveries, j)
      if is_true[j] == 1 { true_pos = true_pos + 1 }
      else { false_pos = false_pos + 1 }
    }
  }
  let total = len(discoveries)
  let fdr = if total > 0 then false_pos / total else 0.0
  print("{method_name} | {total:>17} | {true_pos:>8} | {false_pos:>9} | {fdr:>6.3}")
}

Drawing a Volcano Plot

# Volcano plot: the most important visualization in DE analysis
let log2_fc = fold_changes
let neg_log10_p = p_bh |> map(|p| if p > 0 then -log10(p) else 10)

# Classify genes
let colors = []
for i in 0..len(p_bh) {
  if p_bh[i] < 0.05 and abs(log2_fc[i]) > 1.0 {
    colors = append(colors, "significant")
  } else {
    colors = append(colors, "not_significant")
  }
}

volcano(log2_fc, neg_log10_p, {title: "Differential Expression: Tumor vs Normal", x_label: "log2 Fold Change", y_label: "-log10(FDR-adjusted p-value)", fc_threshold: 1.0, p_threshold: 0.05, highlight: colors})

# Count genes in each quadrant
let sig_up = 0
let sig_down = 0
let not_sig = 0
for i in 0..len(p_bh) {
  if p_bh[i] < 0.05 and log2_fc[i] > 1.0 { sig_up = sig_up + 1 }
  else if p_bh[i] < 0.05 and log2_fc[i] < -1.0 { sig_down = sig_down + 1 }
  else { not_sig = not_sig + 1 }
}
print("Significantly upregulated:   {sig_up}")
print("Significantly downregulated: {sig_down}")
print("Not significant:             {not_sig}")

Visualizing the BH Procedure

set_seed(42)
# Step-by-step BH procedure visualization

# Smaller example for clarity: 100 tests, 10 truly different
let p_vals = []
for i in 1..100 {
  let g1 = rnorm(10, 0, 1)
  let g2 = if i <= 10 {
    rnorm(10, 3, 1)
  } else {
    rnorm(10, 0, 1)
  }
  p_vals = append(p_vals, ttest(g1, g2).p_value)
}

# Sort p-values
let sorted_p = sort(p_vals)
let bh_thresholds = range(1, 101) |> map(|i| i / 100 * 0.05)

# Plot sorted p-values against BH thresholds
scatter(range(1, 101), sorted_p, {title: "BH Procedure: Sorted p-values vs Threshold Line", x_label: "Rank", y_label: "p-value", overlay_lines: [[range(1, 101), bh_thresholds]]})

let bh_adjusted = p_adjust(p_vals, "BH")
let n_sig = bh_adjusted |> filter(|p| p < 0.05) |> len()
print("BH discoveries (FDR < 0.05): {n_sig}")

p-Value Histograms: Diagnostic Tool

set_seed(42)
# A well-behaved p-value histogram tells you a lot
# Uniform = all null; spike at 0 = true signal exists


# Scenario 1: All null (should be uniform)
let null_ps = []
for i in 1..5000 {
  let g1 = rnorm(10, 0, 1)
  let g2 = rnorm(10, 0, 1)
  null_ps = append(null_ps, ttest(g1, g2).p_value)
}

# Scenario 2: 20% true signal (spike near 0 + uniform)
let mixed_ps = []
for i in 1..5000 {
  let g1 = rnorm(10, 0, 1)
  let g2 = if i <= 1000 {
    rnorm(10, 2, 1)
  } else {
    rnorm(10, 0, 1)
  }
  mixed_ps = append(mixed_ps, ttest(g1, g2).p_value)
}

histogram(null_ps, {title: "All Null: Uniform p-value Distribution", x_label: "p-value", bins: 20})

histogram(mixed_ps, {title: "20% True Signal: Spike Near Zero + Uniform Background", x_label: "p-value", bins: 20})

Python:

from statsmodels.stats.multitest import multipletests
import numpy as np

# Simulate p-values (example with 1000 tests)
np.random.seed(42)
p_values = np.random.uniform(0, 1, 1000)
p_values[:100] = np.random.uniform(0, 0.01, 100)  # 100 true signals

# Apply corrections
reject_bonf, padj_bonf, _, _ = multipletests(p_values, method='bonferroni')
reject_holm, padj_holm, _, _ = multipletests(p_values, method='holm')
reject_bh, padj_bh, _, _ = multipletests(p_values, method='fdr_bh')

print(f"Bonferroni: {reject_bonf.sum()} discoveries")
print(f"Holm:       {reject_holm.sum()} discoveries")
print(f"BH (FDR):   {reject_bh.sum()} discoveries")

R:

# Simulate p-values
set.seed(42)
p_values <- c(runif(100, 0, 0.01), runif(900, 0, 1))

# Apply corrections (all built-in)
p_bonf <- p.adjust(p_values, method = "bonferroni")
p_holm <- p.adjust(p_values, method = "holm")
p_bh   <- p.adjust(p_values, method = "BH")
p_by   <- p.adjust(p_values, method = "BY")

cat("Bonferroni:", sum(p_bonf < 0.05), "\n")
cat("Holm:      ", sum(p_holm < 0.05), "\n")
cat("BH (FDR):  ", sum(p_bh < 0.05), "\n")
cat("BY:        ", sum(p_by < 0.05), "\n")

# Volcano plot (using EnhancedVolcano)
library(EnhancedVolcano)
EnhancedVolcano(results, lab = rownames(results),
  x = 'log2FoldChange', y = 'padj',
  pCutoff = 0.05, FCcutoff = 1.0)

Exercises

Exercise 1: Simulate and Count

Simulate 10,000 tests where all genes are null (no true differences). Verify that approximately 500 have p < 0.05, 100 have p < 0.01, and 10 have p < 0.001.

# TODO: Simulate 10,000 null t-tests
# TODO: Count p < 0.05, p < 0.01, p < 0.001
# TODO: Plot the p-value histogram (should be uniform)

Exercise 2: Compare Correction Power

Simulate 5,000 genes (500 truly DE with log2FC = 1.5, 4,500 null). Apply Bonferroni, Holm, and BH. For each, compute: (a) total discoveries, (b) true positives, (c) false positives, (d) actual FDR.

# TODO: Simulate 5,000 genes (500 DE + 4,500 null)
# TODO: Apply all three corrections
# TODO: Build a comparison table
# TODO: Which method gives the best balance of power and error control?

Exercise 3: Build a Volcano Plot

Using the simulation from Exercise 2, create a volcano plot. Color genes as “significant” (BH-adjusted p < 0.05 AND |log2FC| > 1) versus “not significant.”

# TODO: Use fold changes and BH-adjusted p-values from Exercise 2
# TODO: Create a volcano plot with appropriate thresholds
# TODO: Count genes in each quadrant

Exercise 4: p-Value Histogram Diagnostics

Generate p-value histograms for three scenarios: (a) 100% null genes, (b) 10% true DE genes, (c) 50% true DE genes. Describe the characteristic shape of each and explain what you would look for in real data.

# TODO: Scenario A — all null
# TODO: Scenario B — 10% DE
# TODO: Scenario C — 50% DE
# TODO: Create histograms for each
# TODO: Describe the shape and what it tells you

Exercise 5: When Does Bonferroni Make Sense?

You are testing 5 pre-specified candidate genes (not genome-wide). Apply Bonferroni and BH to these 5 p-values: [0.008, 0.023, 0.041, 0.062, 0.110]. How do the results differ? When would you prefer Bonferroni here?

let p_vals = [0.008, 0.023, 0.041, 0.062, 0.110]

# TODO: Apply Bonferroni and BH
# TODO: Which genes are significant under each method?
# TODO: Argue for which method is more appropriate for 5 candidate genes

Key Takeaways

• Testing m hypotheses at alpha = 0.05 yields approximately m x 0.05 false positives — devastating for genomics
• Bonferroni correction (p x m) controls FWER but is often too conservative for large-scale studies
• Holm’s step-down method is always at least as powerful as Bonferroni and should be preferred
• Benjamini-Hochberg (BH) FDR correction is the standard for genomics: it controls the expected proportion of false discoveries rather than the probability of any false discovery
• At FDR = 0.05, you accept that about 5% of discoveries may be false — a practical trade-off for high-throughput biology
• p-value histograms are essential diagnostics: uniform = all null, spike at zero = true signal present
• The volcano plot (-log10 adjusted p vs log2 fold change) is the standard visualization for differential expression: genes in the upper corners are both significant and biologically meaningful
• Always report which correction method you used — “p < 0.05” means very different things with and without adjustment

What’s Next

With the multiple testing crisis solved, we have completed the core toolkit of statistical hypothesis testing. Week 3 shifts to modeling and relationships: tomorrow we explore correlation — how to measure and test whether two continuous variables move together, from gene expression co-regulation to dose-response curves. This transition from “are groups different?” to “how are variables related?” opens the door to regression, prediction, and the modeling approaches that power modern biostatistics.

Day 13: Correlation — Finding Relationships

The Problem

Dr. Sarah Kim is studying breast cancer transcriptomics across 200 tumor samples. She notices that BRCA1 and BARD1 expression seem to rise and fall together — when one is high, the other tends to be high too. Exciting! These genes encode proteins that form a heterodimer critical for DNA repair.

But her collaborator raises a concern: “Both genes are upregulated in rapidly dividing cells. Couldn’t cell proliferation be driving both signals independently? You might be seeing a spurious association.”

Sarah needs to:

1. Quantify how strongly BRCA1 and BARD1 co-vary
2. Determine whether the relationship is statistically significant
3. Control for the confounding effect of proliferation
4. Visualize relationships across an entire panel of DNA repair genes

This is the domain of correlation analysis — one of the most used (and most misused) tools in all of biology.

What Is Correlation?

Correlation measures the strength and direction of the relationship between two variables. It answers: “When one variable goes up, does the other tend to go up, go down, or do nothing?”

The correlation coefficient ranges from -1 to +1:

Key insight: Correlation is symmetric — the correlation between X and Y is the same as between Y and X. It doesn’t imply direction or causation.

Three Types of Correlation

1. Pearson’s r: The Linear Standard

Pearson’s correlation coefficient measures linear association:

$$r = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i - \bar{x})^2 \cdot \sum_{i=1}^{n}(y_i - \bar{y})^2}}$$

Assumptions:

• Both variables are continuous
• Relationship is approximately linear
• Data are roughly normally distributed (for inference)
• No extreme outliers (a single outlier can flip the sign)

Strengths: Most powerful when assumptions are met. Directly related to R² in linear regression.

Weakness: Sensitive to outliers. Misses non-linear relationships entirely.

2. Spearman’s ρ (rho): The Rank-Based Alternative

Spearman’s correlation operates on ranks rather than raw values. It measures monotonic association — whether Y tends to increase (or decrease) as X increases, even if not linearly.

How it works:

1. Rank each variable separately (1, 2, 3, …)
2. Compute Pearson’s r on the ranks

Assumptions:

• Ordinal or continuous data
• Monotonic relationship (doesn’t need to be linear)

Strengths: Robust to outliers. Works with non-normal distributions. Handles log-scale expression data naturally.

Weakness: Less powerful than Pearson when linearity holds.

3. Kendall’s τ (tau): The Most Robust

Kendall’s tau counts concordant and discordant pairs:

$$\tau = \frac{(\text{concordant pairs}) - (\text{discordant pairs})}{\binom{n}{2}}$$

A pair (i, j) is concordant if both xᵢ > xⱼ and yᵢ > yⱼ (or both less). It’s discordant if they disagree.

Strengths: Most robust to outliers. Better for small samples. Has clearer probabilistic interpretation.

Weakness: Computationally slower for large datasets. Values tend to be smaller than Pearson/Spearman for the same data.

Decision Table: Which Correlation to Use?

Common pitfall: Using Pearson on raw RNA-seq counts. These are heavily right-skewed — a few highly expressed genes dominate the correlation. Always log-transform first or use Spearman.

Partial Correlation: Controlling for Confounders

Standard correlation between X and Y can be inflated or deflated by a confounding variable Z that influences both. Partial correlation removes the effect of Z:

$$r_{XY \cdot Z} = \frac{r_{XY} - r_{XZ} \cdot r_{YZ}}{\sqrt{(1 - r_{XZ}^2)(1 - r_{YZ}^2)}}$$

Example: BRCA1 and BARD1 might correlate at r = 0.85. But if cell proliferation (MKI67 expression) drives both, the partial correlation controlling for MKI67 might drop to r = 0.45 — still real, but weaker.

Clinical relevance: In pharmacogenomics, two drug targets may appear correlated simply because both are overexpressed in a particular cancer subtype. Partial correlation controlling for subtype reveals whether the targets have an independent relationship.

Anscombe’s Quartet: Why You Must Visualize

In 1973, Francis Anscombe created four datasets with identical Pearson correlations (r = 0.816), identical means, and identical regression lines — but wildly different patterns:

The lesson: Never trust a correlation coefficient without a scatter plot.

Correlation Matrix and Heatmaps

When studying many variables simultaneously, a correlation matrix shows all pairwise correlations. For p variables, this is a p × p symmetric matrix with 1s on the diagonal.

For genomics, this is invaluable for:

• Identifying co-expression modules
• Detecting batch effects (technical variables cluster together)
• Revealing pathway relationships

A heatmap visualization uses color intensity to represent correlation strength, making patterns immediately visible across dozens or hundreds of variables.

Testing Significance: Is This Correlation Real?

A correlation of r = 0.3 might be noise in 10 samples but highly significant in 1000. The correlation test evaluates:

• H₀: ρ = 0 (no correlation in the population)
• H₁: ρ ≠ 0

The test statistic follows a t-distribution with n-2 degrees of freedom:

$$t = r\sqrt{\frac{n-2}{1-r^2}}$$

Key insight: With large n (common in genomics), even tiny correlations become “significant.” A correlation of r = 0.05 is significant at p < 0.05 when n > 1500. Always report both the coefficient AND the p-value, and judge practical significance by the magnitude of r.

Correlation ≠ Causation

This cannot be overstated. Correlation tells you variables co-vary — nothing more.

Classic examples in biology:

• Ice cream sales correlate with drowning deaths (confounder: hot weather)
• Shoe size correlates with reading ability in children (confounder: age)
• Stork population correlates with birth rate across European countries (confounder: rural vs. urban)

To establish causation, you need:

1. Temporal precedence — cause precedes effect
2. Experimental manipulation — perturb X, observe Y change
3. Elimination of confounders — no third variable explains both
4. Biological mechanism — plausible pathway

Correlation in BioLang

Basic Correlations

set_seed(42)
# Generate tumor expression data for 200 samples
let n = 200

# Simulate BRCA1 and BARD1 with true correlation + noise
let proliferation = rnorm(n, 10, 3)
let brca1 = proliferation * 0.8 + rnorm(n, 5, 2)
let bard1 = proliferation * 0.7 + rnorm(n, 4, 2)

# Pearson correlation — assumes linearity
let r_pearson = cor(brca1, bard1)
print("Pearson r: {r_pearson}")  # ~0.78

# Spearman rank correlation — monotonic, robust
let r_spearman = spearman(brca1, bard1)
print("Spearman ρ: {r_spearman}")  # ~0.76

# Kendall tau — concordant/discordant pairs, most robust
let r_kendall = kendall(brca1, bard1)
print("Kendall τ: {r_kendall}")  # ~0.57 (typically smaller)

Statistical Testing

# Test whether correlation is significantly different from zero
# cor() returns the coefficient; use a t-test transformation for p-value
let r = cor(brca1, bard1)
let t_stat = r * sqrt((n - 2) / (1 - r * r))
print("r = {r}, t = {t_stat}")

# Spearman returns {coefficient, pvalue}
let spearman_result = spearman(brca1, bard1)
print("Spearman ρ = {spearman_result}")

Partial Correlation: Removing Confounders

set_seed(42)
# BRCA1-BARD1 correlation controlling for proliferation (MKI67)
let mki67 = proliferation + rnorm(n, 0, 1)

# Raw correlation
let r_raw = cor(brca1, bard1)
print("Raw Pearson r: {r_raw}")  # ~0.78

# Partial correlation controlling for MKI67
# Compute manually: regress out confounder from both variables
let r_xz = cor(brca1, mki67)
let r_yz = cor(bard1, mki67)
let r_partial = (r_raw - r_xz * r_yz) / sqrt((1 - r_xz * r_xz) * (1 - r_yz * r_yz))
print("Partial r (controlling MKI67): {r_partial}")  # lower — confounder removed

# The difference reveals how much of the BRCA1-BARD1
# association was driven by shared proliferation signal
print("Reduction: {((r_raw - r_partial) / r_raw * 100) |> round(1)}%")

Correlation Matrix and Heatmap

set_seed(42)
# Build expression matrix for 8 DNA repair genes
let base = rnorm(200, 0, 1)

let genes = {
    "BRCA1":  base * 0.8 + rnorm(200, 10, 2),
    "BARD1":  base * 0.7 + rnorm(200, 8, 2),
    "RAD51":  base * 0.6 + rnorm(200, 12, 3),
    "PALB2":  base * 0.5 + rnorm(200, 9, 2),
    "ATM":    rnorm(200, 11, 3),
    "TP53":   base * -0.4 + rnorm(200, 15, 4),
    "MDM2":   base * -0.3 + rnorm(200, 7, 2),
    "GAPDH":  rnorm(200, 20, 1)
}

# Compute pairwise correlations
let gene_names = ["BRCA1", "BARD1", "RAD51", "PALB2", "ATM", "TP53", "MDM2", "GAPDH"]
for i in 0..8 {
    for j in (i+1)..8 {
        let r = cor(genes[gene_names[i]], genes[gene_names[j]])
        print("{gene_names[i]} vs {gene_names[j]}: r = {r |> round(3)}")
    }
}

# Visualize as heatmap — instantly reveals co-expression modules
heatmap(genes, {title: "DNA Repair Gene Co-Expression", color_scale: "RdBu"})

Visualizing with Scatter Plots

# Scatter plot with correlation annotation
let scatter_data = table({"BRCA1": brca1, "BARD1": bard1})
plot(scatter_data, {type: "scatter", x: "BRCA1", y: "BARD1",
    title: "BRCA1 vs BARD1 Co-Expression",
    x_label: "BRCA1 Expression (log2 FPKM)",
    y_label: "BARD1 Expression (log2 FPKM)"})

Demonstrating Anscombe’s Quartet Effect

set_seed(42)
# Two datasets with identical Pearson r but different patterns
let x_linear = rnorm(100, 10, 3)
let y_linear = x_linear * 0.5 + rnorm(100, 0, 2)

let x_curve = rnorm(100, 10, 5)
let y_curve = (x_curve - 10) ** 2 / 10 + rnorm(100, 0, 1)

print("Linear: Pearson r = {cor(x_linear, y_linear)}")
print("Curved: Pearson r = {cor(x_curve, y_curve)}")
print("Linear: Spearman ρ = {spearman(x_linear, y_linear)}")
print("Curved: Spearman ρ = {spearman(x_curve, y_curve)}")

# Spearman catches the monotonic but non-linear pattern
# Always plot your data!

Python:

import numpy as np
from scipy import stats
import seaborn as sns
import matplotlib.pyplot as plt

# Pearson
r, p = stats.pearsonr(brca1, bard1)

# Spearman
rho, p = stats.spearmanr(brca1, bard1)

# Kendall
tau, p = stats.kendalltau(brca1, bard1)

# Partial correlation (using pingouin)
import pingouin as pg
partial = pg.partial_corr(data=df, x='BRCA1', y='BARD1', covar='MKI67')

# Correlation matrix heatmap
corr = df.corr(method='spearman')
sns.heatmap(corr, annot=True, cmap='RdBu_r', center=0)

R:

# Pearson
cor.test(brca1, bard1, method = "pearson")

# Spearman
cor.test(brca1, bard1, method = "spearman")

# Kendall
cor.test(brca1, bard1, method = "kendall")

# Partial correlation (using ppcor)
library(ppcor)
pcor.test(brca1, bard1, mki67)

# Correlation matrix heatmap
library(corrplot)
corrplot(cor(gene_matrix, method = "spearman"),
         method = "color", type = "upper")

Exercises

Exercise 1: Compare Three Methods

Compute Pearson, Spearman, and Kendall correlations for the following gene pairs. Which method gives the most different result, and why?

set_seed(42)
let n = 150

# Gene pair 1: linear relationship with outliers
let gene_a = rnorm(n, 8, 2)
let gene_b = gene_a * 0.6 + rnorm(n, 3, 1)

# Add 5 extreme outliers
# (imagine contaminated samples)

# Compute all three correlations for gene_a vs gene_b
# Which is most affected by outliers?

Exercise 2: Partial Correlation in Drug Response

Three variables are measured across 100 cancer cell lines: drug sensitivity (IC50), target gene expression, and cell doubling time. The target gene and IC50 appear correlated. Is the relationship real, or driven by growth rate?

set_seed(42)
let n = 100
let growth_rate = rnorm(n, 24, 6)

let target_expr = growth_rate * 0.5 + rnorm(n, 10, 3)
let ic50 = growth_rate * -0.3 + rnorm(n, 50, 10)

# 1. Compute raw correlation between target_expr and ic50
# 2. Compute partial correlation controlling for growth_rate
# 3. What fraction of the apparent association was confounded?

Exercise 3: Build and Interpret a Heatmap

Create a correlation heatmap for the following gene expression panel. Identify which genes cluster into co-expression modules.

set_seed(42)
let n = 300

# Simulate 3 biological modules:
# Module 1 (immune): CD8A, GZMB, PRF1, IFNG
# Module 2 (proliferation): MKI67, TOP2A, PCNA
# Module 3 (housekeeping): GAPDH, ACTB

# Create correlated expression within modules
# and weak/no correlation between modules
# Compute pairwise cor() and visualize with heatmap
# Which modules emerge from the clustering?

Exercise 4: Significance vs. Magnitude

Generate datasets with n = 20 and n = 2000. Show that a weak correlation (r ≈ 0.1) is non-significant with small n but highly significant with large n. Argue why the p-value alone is misleading.

set_seed(42)
# Small sample: n = 20, weak correlation
let x_small = rnorm(20, 0, 1)
let y_small = x_small * 0.1 + rnorm(20, 0, 1)

# Large sample: n = 2000, same weak correlation
let x_large = rnorm(2000, 0, 1)
let y_large = x_large * 0.1 + rnorm(2000, 0, 1)

# Compute cor() on both and test significance
# Compare p-values and correlation magnitudes
# What should you report?

Exercise 5: Anscombe Challenge

Create two synthetic gene-pair datasets where Pearson r is nearly identical (~0.7) but scatter plots reveal completely different biology — one linear, one with a threshold effect (flat then rising). Show that Spearman catches the difference.

Key Takeaways

• Pearson measures linear association; Spearman measures monotonic (rank-based); Kendall counts concordant pairs and is most robust
• Correlation ranges from -1 to +1; the sign indicates direction, the magnitude indicates strength
• Always visualize — identical correlations can hide wildly different patterns (Anscombe’s quartet)
• Partial correlation removes the effect of confounders, revealing true associations
• With large genomics datasets, tiny correlations become “significant” — always report the magnitude alongside the p-value
• Correlation is not causation — co-expression does not imply co-regulation or functional relationship
• Spearman is generally the safest default for gene expression data

What’s Next

Now that we can quantify relationships between two variables, Day 14 takes the next step: using one variable to predict another. We’ll build our first linear regression models, learning to fit lines, interpret slopes, and critically evaluate whether our predictions are trustworthy.

Day 14: Linear Regression — Prediction from Data

The Problem

Dr. James Park is a pharmacogenomicist working with the NCI-60 cell line panel — 60 cancer cell lines spanning 9 tissue types. He has gene expression data and drug sensitivity measurements (IC50 values) for each line. His question: Can we predict how sensitive a cell line will be to a new kinase inhibitor based on its expression of the drug’s target gene?

He knows the target gene and IC50 seem correlated (Day 13 confirmed r = -0.72). But correlation just says “they move together.” James needs to go further: given a new cell line with a known expression level, what IC50 should he predict? And how confident should he be?

This is the leap from association to prediction — the domain of linear regression.

What Is Linear Regression?

Linear regression fits a straight line through data points to model the relationship between a predictor (X) and a response (Y):

$$Y = \beta_0 + \beta_1 X + \varepsilon$$

The method of least squares finds the β₀ and β₁ that minimize the sum of squared residuals:

$$\min_{\beta_0, \beta_1} \sum_{i=1}^{n} (y_i - \beta_0 - \beta_1 x_i)^2$$

Key insight: Regression has a clear asymmetry that correlation lacks. X predicts Y — there is a designated predictor and a designated outcome. The regression of Y on X is NOT the same as X on Y.

From Correlation to Regression

Correlation and simple regression are intimately linked:

If r = -0.72, then R² = 0.52, meaning 52% of the variation in IC50 is explained by target expression. The remaining 48% is unexplained — due to other genes, pathway redundancy, or noise.

Interpreting Regression Output

A regression produces a table of coefficients:

Reading this:

• Intercept: When expression = 0, predicted IC50 is 85.3 μM (often not biologically meaningful)
• Slope: Each 1-unit increase in expression decreases IC50 by 4.7 μM (more expression → more sensitive)
• p-value: The slope is significantly different from zero — expression truly predicts sensitivity
• R²: How well the line fits overall

Common pitfall: The intercept often represents an extrapolation beyond the data range. If expression ranges from 5-15, interpreting “IC50 when expression = 0” is meaningless. Focus on the slope.

Prediction: Point Estimates and Intervals

Once we have a fitted model, we can predict Y for new X values. But predictions come with two types of uncertainty:

Confidence Interval (for the mean response): “We’re 95% confident the average IC50 for all cell lines with expression = 10 lies in this range.”

Prediction Interval (for a single new observation): “We’re 95% confident the IC50 for one specific new cell line with expression = 10 lies in this range.”

Prediction intervals are always wider than confidence intervals because they include individual-level noise (ε).

Clinical relevance: In precision medicine, you’re usually predicting for one patient (prediction interval), not the population average (confidence interval). The prediction interval honestly reflects your uncertainty.

Residual Analysis: Checking Your Model

A regression model makes assumptions. Residuals (observed - predicted) reveal violations:

The four assumptions of linear regression:

1. Linearity: Y is a linear function of X (check scatter plot)
2. Independence: Observations are independent (study design)
3. Normality: Residuals are normally distributed (Q-Q plot)
4. Homoscedasticity: Constant error variance (residuals vs. fitted)

Common pitfall: Many biological relationships are not linear on the original scale but become linear after log-transformation. Always consider transforming IC50 or expression values before fitting.

When Regression Goes Wrong

Extrapolation

Predicting beyond the range of your data is dangerous. If expression ranges from 5-15 in your training data, predicting IC50 for expression = 25 assumes the linear trend continues — it may not.

Confounding

A significant regression doesn’t prove causation. The apparent effect of expression on IC50 could be mediated by a third variable (e.g., tissue type).

Non-linearity

If the true relationship is curved (threshold effect, saturation), a line fits poorly. Residual plots expose this.

Influential Points

A single extreme observation can dramatically change the fitted line. Leverage and Cook’s distance help identify such points.

Linear Regression in BioLang

Simple Linear Regression

set_seed(42)
# NCI-60 pharmacogenomics: predict IC50 from target expression
let n = 60

# Simulate expression and drug sensitivity
let expression = rnorm(n, 10, 3)
let ic50 = 85 - 4.5 * expression + rnorm(n, 0, 8)

# Fit simple linear regression
let model = lm(ic50, expression)

# Print model summary
print("=== Linear Regression Summary ===")
print("Intercept: {model.intercept}")
print("Slope: {model.slope}")
print("R²: {model.r_squared}")
print("Adjusted R²: {model.adj_r_squared}")
print("F-statistic p-value: {model.p_value}")

Interpreting Coefficients

# Detailed coefficient table
print("=== Coefficients ===")
print("  Intercept: {model.intercept}")
print("  Expression β: {model.slope}")
print("  p-value: {model.p_value}")

# Interpretation
let slope = model.slope
print("\nFor each 1-unit increase in expression,")
print("IC50 decreases by {slope |> abs |> round(2)} μM")
print("R² = {model.r_squared |> round(3)}: expression explains")
print("{(model.r_squared * 100) |> round(1)}% of IC50 variation")

Prediction with Intervals

# Predict IC50 for new cell lines
let new_expression = [7.0, 10.0, 13.0]

# Point predictions using model coefficients
print("=== Predictions ===")
for i in 0..3 {
    let predicted = model.slope * new_expression[i] + model.intercept
    print("Expression = {new_expression[i]}: Predicted IC50 = {predicted |> round(1)} μM")
}

Residual Analysis

# Compute residuals and fitted values
let fitted = expression |> map(|x| model.slope * x + model.intercept)
let resid = []
for i in 0..n {
    resid = resid + [ic50[i] - fitted[i]]
}

# 1. Residuals vs Fitted — check for patterns
let resid_table = table({"Fitted": fitted, "Residual": resid})
plot(resid_table, {type: "scatter", x: "Fitted", y: "Residual",
    title: "Residuals vs Fitted"})

# 2. Check residual distribution
print("Residual summary:")
print(summary(resid))

Visualization: Scatter with Regression Line

# Scatter plot with regression line
let plot_data = table({"Expression": expression, "IC50": ic50})
plot(plot_data, {type: "scatter", x: "Expression", y: "IC50",
    title: "Drug Sensitivity vs Target Expression (NCI-60)",
    x_label: "Target Gene Expression (log2 FPKM)",
    y_label: "Drug IC50 (μM)"})

Demonstrating Problems: Extrapolation

set_seed(42)
# Show danger of extrapolation
let x = rnorm(50, 10, 3)
let y = 100 - 3 * x + 0.2 * x ** 2 + rnorm(50, 0, 3)

# Fit linear model in observed range
let model_extrap = lm(y, x)
print("R² in training range: {model_extrap.r_squared |> round(3)}")

# Predict within range — reasonable
let pred_in = model_extrap.slope * 10.0 + model_extrap.intercept
print("Prediction at x=10 (in range): {pred_in |> round(1)}")

# Predict outside range — dangerous!
let pred_out = model_extrap.slope * 25.0 + model_extrap.intercept
print("Prediction at x=25 (extrapolation): {pred_out |> round(1)}")
print("True value at x=25 would be: {(100 - 3*25 + 0.2*25**2) |> round(1)}")
print("Extrapolation error demonstrates the danger!")

Working with Log-Transformed Data

set_seed(42)
# Many biological relationships are linear on the log scale
let dose = rnorm(40, 50, 25) |> map(|d| max(d, 0.1))
let response = 50 / (1 + (dose / 10) ** 1.5) + rnorm(40, 0, 3)

# Linear model on raw scale — poor fit
let model_raw = lm(response, dose)
print("R² (raw scale): {model_raw.r_squared |> round(3)}")

# Log-transform dose — much better
let log_dose = dose |> map(|d| log2(d))
let model_log = lm(response, log_dose)
print("R² (log scale): {model_log.r_squared |> round(3)}")

# Compare residual plots to see the difference

Python:

import numpy as np
from scipy import stats
import statsmodels.api as sm

# Simple linear regression
X = sm.add_constant(expression)  # adds intercept
model = sm.OLS(ic50, X).fit()
print(model.summary())

# Prediction with intervals
predictions = model.get_prediction(sm.add_constant([7, 10, 13]))
pred_summary = predictions.summary_frame(alpha=0.05)
# obs_ci_lower/obs_ci_upper = prediction interval
# mean_ci_lower/mean_ci_upper = confidence interval

# Residual analysis
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].scatter(model.fittedvalues, model.resid)
axes[0].axhline(0, color='red')
stats.probplot(model.resid, plot=axes[1])

R:

# Simple linear regression
model <- lm(ic50 ~ expression)
summary(model)

# Prediction with intervals
new_data <- data.frame(expression = c(7, 10, 13))
predict(model, new_data, interval = "confidence")
predict(model, new_data, interval = "prediction")

# Residual diagnostics — 4 diagnostic plots
par(mfrow = c(2, 2))
plot(model)

Exercises

Exercise 1: Fit and Interpret

A study measures tumor mutation burden (TMB) and neoantigen count across 80 melanoma samples. Fit a linear regression and answer: How many additional neoantigens does each additional mutation predict?

set_seed(42)
let n = 80
let tmb = rnorm(n, 200, 80)
let neoantigens = 5 + 0.15 * tmb + rnorm(n, 0, 12)

# 1. Fit lm(neoantigens, tmb)
# 2. What is the slope? Interpret it biologically
# 3. What is model.r_squared? Is TMB a good predictor of neoantigen count?
# 4. Predict neoantigen count for TMB = 100, 200, 400

Exercise 2: Residual Diagnostics

Fit a linear model to the dose-response data below and use residual plots to determine whether the linear assumption holds. If not, what transformation improves the fit?

set_seed(42)
let dose = rnorm(60, 25, 12) |> map(|d| max(d, 1))
let effect = 20 * log2(dose) + rnorm(60, 0, 5)

# 1. Fit lm(effect, dose)
# 2. Create residuals vs fitted plot — what pattern do you see?
# 3. Try log-transforming dose and re-fitting
# 4. Compare R² values and residual patterns

Exercise 3: Confidence vs. Prediction Intervals

Demonstrate why prediction intervals are always wider than confidence intervals. Generate data, fit a model, and plot both intervals across the predictor range.

set_seed(42)
let x = rnorm(50, 10, 5)
let y = 10 + 2 * x + rnorm(50, 0, 4)

# 1. Fit lm(y, x)
# 2. Generate predictions for x = 0, 2, 4, ..., 20
# 3. Use model.slope * xi + model.intercept for each
# 4. Compare point estimates at different x values

Exercise 4: The Outlier Effect

Add a single influential outlier to well-behaved data and show how it changes the slope, R², and predictions. Then remove it and compare.

set_seed(42)
let x_clean = rnorm(49, 10, 2)
let y_clean = 5 + 1.5 * x_clean + rnorm(49, 0, 2)

# Add one extreme outlier at x=10, y=50 (should be ~20)
# 1. Fit model with and without the outlier
# 2. Compare slopes and R² values
# 3. How much does prediction at x=12 change?

Key Takeaways

• Linear regression models Y = β₀ + β₁X + ε, using least squares to find the best-fit line
• R² tells you the proportion of variance explained; the slope tells you the rate of change
• Prediction intervals (for individuals) are always wider than confidence intervals (for the mean) — use the right one for your question
• Residual analysis is mandatory: check linearity, normality, and constant variance before trusting results
• Log-transforming variables often linearizes biological relationships (dose-response, expression-phenotype)
• Beware extrapolation — the linear trend may not continue beyond your data range
• A significant regression does not prove causation — confounders may drive the relationship
• Always visualize your data with a scatter plot before and after fitting

What’s Next

One predictor is rarely enough. Day 15 introduces multiple regression, where we predict outcomes from many variables simultaneously — and face new challenges like multicollinearity, model selection, and regularization.

Day 15: Multiple Regression and Model Selection

The Problem

Dr. Maria Chen is a clinical researcher studying pancreatic cancer. She has tumor samples from 120 patients, each profiled with 10 biomarkers: CA19-9, CEA, MKI67, TP53 status, tumor size, age, albumin, CRP, neutrophil-lymphocyte ratio (NLR), and platelet count. She wants to predict tumor stage (a continuous composite score from 1.0 to 4.0) from these biomarkers.

But there’s a problem. CA19-9 and CEA are highly correlated (r = 0.88) — they measure overlapping biology. Including both inflates standard errors and makes coefficients uninterpretable. And with 10 potential predictors, how does she find the best subset without overfitting?

She needs multiple regression with careful model selection.

What Is Multiple Regression?

Multiple regression extends simple regression to multiple predictors:

$$Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p + \varepsilon$$

Each coefficient βⱼ represents the effect of Xⱼ holding all other predictors constant. This is fundamentally different from running p separate simple regressions.

Key insight: In simple regression, tumor size might predict stage with β = 0.8. In multiple regression controlling for CA19-9, the tumor size coefficient might drop to β = 0.3 — because CA19-9 captures much of the same information.

Multicollinearity: When Predictors Are Too Similar

Multicollinearity occurs when predictors are highly correlated with each other. It doesn’t bias predictions, but it wreaks havoc on coefficient interpretation:

Detecting Multicollinearity: VIF

The Variance Inflation Factor quantifies how much each predictor is explained by the others:

$$VIF_j = \frac{1}{1 - R_j^2}$$

where R²ⱼ is the R² from regressing Xⱼ on all other predictors.

Common pitfall: VIF > 10 doesn’t always mean “drop the variable.” In genomics, biologically meaningful predictors may be correlated. Consider combining them (e.g., principal component) or using regularization instead.

Model Selection: Finding the Best Model

With p predictors, there are 2ᵖ possible models. For p = 10, that’s 1024 models. We need principled ways to choose.

Information Criteria

Where L = likelihood, k = number of parameters, n = sample size.

AIC tends to select slightly larger models (better prediction). BIC tends to select smaller models (better interpretation). When they disagree, consider your goal.

Stepwise Regression

Automated search through predictor combinations:

Common pitfall: Stepwise regression is exploratory, not confirmatory. The selected model may not replicate. Always validate on held-out data or use cross-validation.

Regularized Regression: Handling Many Predictors

When you have many predictors (especially in genomics where p >> n), ordinary least squares fails. Regularization adds a penalty term:

Ridge Regression (L2)

$$\min \sum(y_i - \hat{y}_i)^2 + \lambda \sum \beta_j^2$$

• Shrinks coefficients toward zero but never to exactly zero
• Handles multicollinearity gracefully
• Good when all predictors might be relevant

Lasso Regression (L1)

$$\min \sum(y_i - \hat{y}_i)^2 + \lambda \sum |\beta_j|$$

• Can shrink coefficients to exactly zero (automatic variable selection)
• Produces sparse models (easy to interpret)
• Preferred when you suspect only a few predictors matter

Elastic Net

$$\min \sum(y_i - \hat{y}_i)^2 + \lambda_1 \sum |\beta_j| + \lambda_2 \sum \beta_j^2$$

• Combines L1 and L2 penalties
• Good when predictors are correlated groups (keeps all or drops all from a group)
• Controlled by mixing parameter α (0 = ridge, 1 = lasso)

Polynomial Regression

When the relationship is curved but you want to stay in the regression framework:

$$Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \varepsilon$$

Useful for dose-response curves, growth curves, and non-linear biomarker relationships. But beware overfitting with high-degree polynomials.

Clinical relevance: The relationship between BMI and mortality is U-shaped — both low and high BMI increase risk. A linear model misses this entirely; a quadratic model captures it.

Multiple Regression in BioLang

Building a Multiple Regression Model

set_seed(42)
# Pancreatic cancer: predict tumor stage from biomarkers
let n = 120

# Simulate correlated biomarkers
let age = rnorm(n, 65, 10)
let tumor_size = rnorm(n, 3.5, 1.2)
let ca19_9 = tumor_size * 50 + rnorm(n, 100, 40)
let cea = ca19_9 * 0.3 + rnorm(n, 5, 3)  # correlated with CA19-9
let mki67 = rnorm(n, 30, 15)
let albumin = rnorm(n, 3.5, 0.5)
let crp = rnorm(n, 15, 10)
let nlr = rnorm(n, 4, 2)

# True model: stage depends on tumor_size, ca19_9, mki67, age
let stage = 1.0 + 0.4 * tumor_size + 0.003 * ca19_9 + 0.01 * mki67
    + 0.01 * age - 0.3 * albumin
    + rnorm(n, 0, 0.3)

# Fit multiple regression with all predictors
let data = table({
    "stage": stage, "age": age, "tumor_size": tumor_size,
    "ca19_9": ca19_9, "cea": cea, "mki67": mki67,
    "albumin": albumin, "crp": crp, "nlr": nlr
})
let model_full = lm("stage ~ age + tumor_size + ca19_9 + cea + mki67 + albumin + crp + nlr", data)

print("=== Full Model ===")
print("R²: {model_full.r_squared |> round(3)}")
print("Adjusted R²: {model_full.adj_r_squared |> round(3)}")

Checking Multicollinearity

# Check multicollinearity via pairwise correlations
let pred_names = ["age", "tumor_size", "ca19_9", "cea", "mki67",
                  "albumin", "crp", "nlr"]
let predictors = [age, tumor_size, ca19_9, cea, mki67, albumin, crp, nlr]

print("=== Pairwise Correlations (VIF proxy) ===")
for i in 0..8 {
    for j in (i+1)..8 {
        let r = cor(predictors[i], predictors[j])
        if abs(r) > 0.7 {
            print("  {pred_names[i]} vs {pred_names[j]}: r = {r |> round(3)} *** HIGH")
        }
    }
}

# CA19-9 and CEA likely show high correlation

Stepwise Model Selection

# Manual model comparison: fit reduced models and compare R²
# Drop CEA (collinear with CA19-9) and noise variables (CRP, NLR)
let data_reduced = table({
    "stage": stage, "age": age, "tumor_size": tumor_size,
    "ca19_9": ca19_9, "mki67": mki67, "albumin": albumin
})
let model_reduced = lm("stage ~ age + tumor_size + ca19_9 + mki67 + albumin", data_reduced)

print("=== Model Comparison ===")
print("Full model R²:    {model_full.r_squared |> round(3)}")
print("Full model Adj R²:    {model_full.adj_r_squared |> round(3)}")
print("Reduced model R²: {model_reduced.r_squared |> round(3)}")
print("Reduced model Adj R²: {model_reduced.adj_r_squared |> round(3)}")
print("If Adj R² is similar, the simpler model is preferred")

Regularized Regression

# Regularized regression concepts
# Note: Ridge, Lasso, and Elastic Net are advanced methods
# typically used when p >> n (many predictors, few samples).
# BioLang provides lm() for standard regression; for regularization,
# use Python (scikit-learn) or R (glmnet) as shown below.

# Demonstrate the concept: compare full vs sparse models
# A "lasso-like" approach: fit models dropping one predictor at a time
# and see which predictors contribute the least
let predictors_list = ["age", "tumor_size", "ca19_9", "cea", "mki67", "albumin", "crp", "nlr"]

print("=== Variable Importance (drop-one analysis) ===")
let full_r2 = model_full.r_squared
print("Full model R²: {full_r2 |> round(4)}")

# Compare by dropping noise predictors
let model_no_crp = lm("stage ~ age + tumor_size + ca19_9 + cea + mki67 + albumin + nlr", data)
let model_no_nlr = lm("stage ~ age + tumor_size + ca19_9 + cea + mki67 + albumin + crp", data)
print("Without CRP: R² = {model_no_crp.r_squared |> round(4)}")
print("Without NLR: R² = {model_no_nlr.r_squared |> round(4)}")
print("Predictors with minimal R² drop are candidates for removal")

Polynomial Regression

set_seed(42)
# Non-linear biomarker relationship
let bmi = rnorm(100, 29, 5)
let risk = 0.5 + 0.1 * (bmi - 25) ** 2 + rnorm(100, 0, 2)

# Linear fit — misses the U-shape
let model_linear = lm(risk, bmi)
print("Linear R²: {model_linear.r_squared |> round(3)}")

# Polynomial fit — add bmi² term
let bmi_sq = bmi |> map(|x| x ** 2)
let poly_data = table({"risk": risk, "bmi": bmi, "bmi_sq": bmi_sq})
let model_poly = lm("risk ~ bmi + bmi_sq", poly_data)
print("Quadratic R²: {model_poly.r_squared |> round(3)}")

# Visualize the improvement
let plot_data = table({"BMI": bmi, "Risk": risk})
plot(plot_data, {type: "scatter", x: "BMI", y: "Risk",
    title: "BMI vs Risk: Linear vs Quadratic Fit"})

Predicted vs. Actual Plot

# The ultimate model validation plot
# Compute predicted values from the reduced model
let predicted_stage = model_reduced.fitted

let pred_actual = table({"Actual": stage, "Predicted": predicted_stage})
plot(pred_actual, {type: "scatter", x: "Actual", y: "Predicted",
    title: "Predicted vs Actual (Reduced Model)"})

# Points along the diagonal = good predictions
# Systematic deviations = model problems
let r_pred = cor(stage, predicted_stage)
print("Correlation (actual vs predicted): {r_pred |> round(3)}")

Python:

import statsmodels.api as sm
from sklearn.linear_model import Lasso, Ridge, ElasticNet
from sklearn.preprocessing import PolynomialFeatures
from statsmodels.stats.outliers_influence import variance_inflation_factor

# Multiple regression
X = sm.add_constant(df[['age', 'tumor_size', 'ca19_9', 'cea', 'mki67']])
model = sm.OLS(stage, X).fit()
print(model.summary())

# VIF
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]

# Stepwise (manual — no built-in in statsmodels)
# Use mlxtend: from mlxtend.feature_selection import SequentialFeatureSelector

# Lasso
from sklearn.linear_model import LassoCV
lasso = LassoCV(cv=5).fit(X_scaled, y)
print(f"Non-zero: {(lasso.coef_ != 0).sum()}")

R:

# Multiple regression
model <- lm(stage ~ age + tumor_size + ca19_9 + cea + mki67 +
            albumin + crp + nlr)
summary(model)

# VIF
library(car)
vif(model)

# Stepwise
step_model <- step(model, direction = "both")

# Lasso
library(glmnet)
cv_lasso <- cv.glmnet(X_matrix, stage, alpha = 1)
coef(cv_lasso, s = "lambda.min")

# Ridge
cv_ridge <- cv.glmnet(X_matrix, stage, alpha = 0)

# Elastic net
cv_enet <- cv.glmnet(X_matrix, stage, alpha = 0.5)

Exercises

Exercise 1: Build and Compare Models

Given 150 samples with 6 predictors, build a full model, then use stepwise selection to find a reduced model. Compare using AIC, BIC, and adjusted R².

set_seed(42)
let n = 150
let x1 = rnorm(n, 10, 3)
let x2 = rnorm(n, 5, 2)
let x3 = x1 * 0.5 + rnorm(n, 0, 1)  # correlated with x1
let x4 = rnorm(n, 20, 5)
let x5 = rnorm(n, 0, 1)  # noise
let x6 = rnorm(n, 0, 1)  # noise

let y = 5 + 2 * x1 + 1.5 * x2 + 0.5 * x4 + rnorm(n, 0, 3)

# 1. Fit full model with all 6 predictors using lm()
# 2. Check pairwise cor() — which predictors are collinear?
# 3. Compare full vs reduced models by adj_r_squared
# 4. Which predictors matter? Do they match the true model?

Exercise 2: Ridge vs. Lasso

Compare ridge and lasso on a dataset where only 3 of 8 predictors truly matter. Which method correctly identifies the true predictors?

set_seed(42)
let n = 100

# 8 predictors, only 3 are real
# Fit lm() with all 8, then with only the true 3
# Compare R² — does the reduced model perform similarly?
# Note: for true ridge/lasso, use Python (scikit-learn) or R (glmnet)

Exercise 3: Polynomial vs. Linear

Fit linear, quadratic, and cubic models to a dose-response curve. Use AIC to select the best model, and show that adding unnecessary complexity (cubic) hurts.

set_seed(42)
let dose = rnorm(80, 25, 12) |> map(|d| max(d, 1))
let response = 10 + 5 * log2(dose) + rnorm(80, 0, 3)

# 1. Fit lm() with dose, dose + dose², dose + dose² + dose³
# 2. Compare R² and adj_r_squared for each
# 3. Which degree gives the best balance of fit and simplicity?

Exercise 4: Predicted vs. Actual

Build a multiple regression model for gene expression prediction. Create a predicted vs. actual scatter plot. Assess where the model succeeds and fails.

set_seed(42)
let n = 200

# Simulate: predict gene expression from 4 TF binding signals
# Build model using lm(), generate predicted vs actual plot
# Calculate mean absolute error
# Identify the 5 worst predictions — what makes them outliers?

Key Takeaways

• Multiple regression estimates the effect of each predictor controlling for all others — fundamentally different from separate simple regressions
• Multicollinearity (VIF > 10) inflates standard errors and makes coefficients uninterpretable — detect it before interpreting
• Model selection balances fit and complexity: AIC favors prediction, BIC favors parsimony, adjusted R² penalizes added terms
• Stepwise regression is useful for exploration but should be validated on held-out data
• Lasso performs automatic variable selection (zeros out irrelevant predictors); Ridge shrinks all coefficients but keeps them; Elastic net combines both
• Polynomial regression captures non-linear relationships within the regression framework but risks overfitting
• Always check residuals, predicted vs. actual plots, and VIF before trusting your model
• With genomics data (p >> n), regularization is not optional — it’s essential

What’s Next

What if your outcome isn’t continuous but binary — responder vs. non-responder, alive vs. dead, mutant vs. wild-type? Day 16 introduces logistic regression, where we predict probabilities of categorical outcomes using ROC curves, odds ratios, and the powerful GLM framework.

Day 16: Logistic Regression — Binary Outcomes

The Problem

Dr. Priya Sharma is an immuno-oncologist analyzing data from 180 melanoma patients who received anti-PD-1 immunotherapy. For each patient, she has three biomarkers — tumor mutational burden (TMB), PD-L1 expression, and microsatellite instability (MSI) status — and one outcome: response (tumor shrank ≥30%) or non-response.

Her first instinct is to use linear regression, predicting response (coded 1/0) from the biomarkers. But the predictions come out as 1.3 for one patient and -0.2 for another. Probabilities can’t be greater than 1 or less than 0.

She needs a method designed for binary outcomes — logistic regression.

Why Linear Regression Fails for Binary Outcomes

When Y is binary (0 or 1), linear regression has fundamental problems:

Key insight: We need a function that maps any real number to the range [0, 1]. The logistic (sigmoid) function does exactly this.

The Logistic Function

The logistic regression model predicts the probability of the outcome being 1:

$$P(Y=1|X) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p)}}$$

The sigmoid function transforms the linear predictor (which ranges from -∞ to +∞) into a probability (which ranges from 0 to 1):

The curve is steepest at P = 0.5 (the decision boundary) and flattens at the extremes — exactly how biological responses behave.

Interpreting Coefficients: Log-Odds and Odds Ratios

Logistic regression coefficients are on the log-odds scale:

$$\ln\left(\frac{P}{1-P}\right) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots$$

This is less intuitive than linear regression. The key transformation:

$$\text{Odds Ratio} = e^{\beta}$$

Key insight: An odds ratio of 2.0 means the odds of response double for each 1-unit increase in the predictor. It does NOT mean the probability doubles — that depends on the baseline probability.

Odds vs. Probability

Common pitfall: When the outcome is common (prevalence > 10%), odds ratios substantially overestimate relative risks. An OR of 3.0 for a 50% baseline event means the probability goes from 50% to 75% — a relative risk of only 1.5. Report both OR and absolute probabilities.

ROC Curves and AUC

The model outputs probabilities. To make yes/no predictions, you need a threshold — above it, predict “responder.” But which threshold?

A Receiver Operating Characteristic (ROC) curve plots sensitivity vs. (1 - specificity) across all possible thresholds:

The Area Under the Curve (AUC) summarizes overall discrimination:

Clinical relevance: In cancer diagnostics, the threshold choice depends on context. Screening tests favor high sensitivity (don’t miss cancers), while confirmatory tests favor high specificity (don’t cause unnecessary biopsies).

The GLM Framework

Logistic regression is a special case of the Generalized Linear Model (GLM) — a flexible framework for non-normal outcomes:

The GLM framework unifies these under one interface, letting you choose the appropriate model for your data type.

Clinical relevance: RNA-seq read counts follow a negative binomial distribution. Tools like DESeq2 use GLMs with a negative binomial family — the same framework as logistic regression, just with a different distribution assumption.

Logistic Regression in BioLang

Fitting a Logistic Model

set_seed(42)
# Immunotherapy response prediction
let n = 180

# Simulate predictors
let tmb = rnorm(n, 10, 5) |> map(|x| max(x, 0))
let pdl1 = rnorm(n, 30, 20) |> map(|x| max(min(x, 100), 0))
let msi = rnorm(n, 0, 1) |> map(|x| if x > 1.0 { 1 } else { 0 })  # ~15% MSI-high

# True response probability (logistic model)
let log_odds = -3 + 0.15 * tmb + 0.02 * pdl1 + 1.5 * msi
let prob = log_odds |> map(|x| 1.0 / (1.0 + exp(-x)))
let response = prob |> map(|p| if rnorm(1, 0, 1)[0] < p { 1 } else { 0 })

print("Response rate: {(response |> sum) / n * 100 |> round(1)}%")

# Fit logistic regression
let glm_data = table({"response": response, "TMB": tmb, "PDL1": pdl1, "MSI": msi})
let model = glm("response ~ TMB + PDL1 + MSI", glm_data, "binomial")

print("=== Logistic Regression Results ===")
print("Intercept: {model.intercept |> round(3)}")

Interpreting Odds Ratios

# Convert coefficients to odds ratios
let coef_names = ["TMB", "PD-L1", "MSI"]
let coefficients = model.coefficients

print("=== Odds Ratios ===")
for i in 0..3 {
    let beta = coefficients[i]
    let or_val = exp(beta)
    print("{coef_names[i]}: β = {beta |> round(3)}, OR = {or_val |> round(2)}")
}

# Interpretation:
# TMB OR = 1.16 means each additional mut/Mb increases odds of response by 16%
# MSI OR = 4.5 means MSI-high patients have 4.5x the odds of responding

ROC Curve and AUC

# Compute predicted probabilities from model
let pred_prob = []
for i in 0..n {
    let lp = model.intercept + model.coefficients[0] * tmb[i]
        + model.coefficients[1] * pdl1[i] + model.coefficients[2] * msi[i]
    pred_prob = pred_prob + [1.0 / (1.0 + exp(-lp))]
}

# ROC curve
let roc_data = table({"actual": response, "predicted": pred_prob})
roc_curve(roc_data)

# Compute AUC
let auc_val = model.auc
print("AUC: {auc_val |> round(3)}")

# Interpretation
if auc_val >= 0.80 {
    print("Good discrimination")
} else if auc_val >= 0.70 {
    print("Acceptable discrimination")
} else {
    print("Poor discrimination — model needs additional predictors")
}

Finding the Optimal Threshold

# Sensitivity/specificity at different thresholds
let thresholds = [0.2, 0.3, 0.4, 0.5, 0.6, 0.7]

print("=== Threshold Analysis ===")
print("Threshold  Sensitivity  Specificity  PPV      NPV")

for t in thresholds {
    let predicted_class = pred_prob |> map(|p| if p >= t { 1 } else { 0 })
    let tp = 0
    let fp = 0
    let tn = 0
    let fn = 0

    for i in 0..n {
        if predicted_class[i] == 1 && response[i] == 1 { tp = tp + 1 }
        if predicted_class[i] == 1 && response[i] == 0 { fp = fp + 1 }
        if predicted_class[i] == 0 && response[i] == 0 { tn = tn + 1 }
        if predicted_class[i] == 0 && response[i] == 1 { fn = fn + 1 }
    }

    let sens = tp / (tp + fn)
    let spec = tn / (tn + fp)
    let ppv = if tp + fp > 0 { tp / (tp + fp) } else { 0 }
    let npv = if tn + fn > 0 { tn / (tn + fn) } else { 0 }

    print("{t}        {sens |> round(3)}       {spec |> round(3)}      {ppv |> round(3)}    {npv |> round(3)}")
}

# Youden's J: optimal balance of sensitivity and specificity

Using the GLM Framework

# Logistic regression via the general GLM interface
let model_glm = glm("response ~ TMB + PDL1 + MSI", glm_data, "binomial")

# Same results, but now you can swap families:

# Poisson regression for count data (e.g., number of mutations)
# let count_data = table({"mutations": mutation_count, "exposure": exposure, "age": age})
# let count_model = glm("mutations ~ exposure + age", count_data, "poisson")

Visualizing Predicted Probabilities

# Box plot: predicted probabilities by actual outcome
let resp_probs = []
let nonresp_probs = []

for i in 0..n {
    if response[i] == 1 {
        resp_probs = resp_probs + [pred_prob[i]]
    } else {
        nonresp_probs = nonresp_probs + [pred_prob[i]]
    }
}

let bp_table = table({"Non-Responders": nonresp_probs, "Responders": resp_probs})
boxplot(bp_table, {title: "Predicted Probabilities by Actual Outcome"})

# Good model: minimal overlap between the two boxes

Effect of Individual Predictors

# Show how each predictor shifts the probability curve
# Fix other predictors at their means
let tmb_range = [0, 5, 10, 15, 20, 25, 30]
let pdl1_mean = 30
let msi_0 = 0

print("=== TMB Effect on Response Probability ===")
print("(PD-L1 = {pdl1_mean}, MSI = negative)")

for t in tmb_range {
    let lp = model.intercept + model.coefficients[0] * t
        + model.coefficients[1] * pdl1_mean + model.coefficients[2] * msi_0
    let p = 1.0 / (1.0 + exp(-lp))
    print("  TMB = {t}: P(response) = {(p * 100) |> round(1)}%")
}

Python:

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_curve, auc, classification_report
import statsmodels.api as sm

# Statsmodels (full inference)
X = sm.add_constant(df[['TMB', 'PDL1', 'MSI']])
model = sm.Logit(response, X).fit()
print(model.summary())

# Odds ratios
print(np.exp(model.params))
print(np.exp(model.conf_int()))

# ROC curve
fpr, tpr, thresholds = roc_curve(response, model.predict(X))
roc_auc = auc(fpr, tpr)
plt.plot(fpr, tpr, label=f'AUC = {roc_auc:.2f}')

# Scikit-learn (prediction-focused)
clf = LogisticRegression().fit(X_train, y_train)
y_prob = clf.predict_proba(X_test)[:, 1]

R:

# Logistic regression
model <- glm(response ~ TMB + PDL1 + MSI, family = binomial)
summary(model)

# Odds ratios with CI
exp(cbind(OR = coef(model), confint(model)))

# ROC curve
library(pROC)
roc_obj <- roc(response, fitted(model))
plot(roc_obj, print.auc = TRUE)

# Optimal threshold (Youden's J)
coords(roc_obj, "best", best.method = "youden")

Exercises

Exercise 1: Build and Interpret a Logistic Model

Predict cancer diagnosis (1 = cancer, 0 = benign) from three biomarkers. Interpret each odds ratio in clinical terms.

set_seed(42)
let n = 200

let biomarker_a = rnorm(n, 50, 15)
let biomarker_b = rnorm(n, 10, 4)
let age = rnorm(n, 55, 12)

let log_odds = -5 + 0.04 * biomarker_a + 0.2 * biomarker_b + 0.03 * age
let prob = log_odds |> map(|x| 1.0 / (1.0 + exp(-x)))
let cancer = prob |> map(|p| if rnorm(1, 0, 1)[0] < p { 1 } else { 0 })

# 1. Fit glm("cancer ~ biomarker_a + biomarker_b + age", data, "binomial")
# 2. Compute and interpret odds ratios: exp(coefficient) for each predictor
# 3. Which biomarker has the strongest effect?
# 4. What does the OR for age mean clinically?

Exercise 2: ROC Analysis

Build a logistic model and evaluate it with an ROC curve. Find the threshold that maximizes Youden’s J index (sensitivity + specificity - 1).

# Use the model from Exercise 1
# 1. Generate predicted probabilities from model coefficients
# 2. Plot ROC curve with roc_curve(table)
# 3. Compute sensitivity and specificity at thresholds 0.3, 0.5, 0.7
# 4. Which threshold maximizes Youden's J?
# 5. If this is a screening test, would you prefer a different threshold? Why?

Exercise 3: Comparing Two Models

Build two logistic models: one with TMB alone, another with TMB + PD-L1 + MSI. Compare their AUC values. Does adding predictors improve discrimination?

set_seed(42)
let n = 150

# Simulate data where TMB is a moderate predictor and MSI adds substantial value
# 1. Fit model_simple = glm("response ~ TMB", data, "binomial")
# 2. Fit model_full = glm("response ~ TMB + PDL1 + MSI", data, "binomial")
# 3. Compare AUC values
# 4. Plot both ROC curves
# 5. Is the improvement worth the added complexity?

Exercise 4: The Separation Problem

What happens when a predictor perfectly separates outcomes? Simulate MSI status that perfectly predicts response and observe what logistic regression does.

set_seed(42)
let n = 100
let msi = rnorm(n, 0, 1) |> map(|x| if x > 0.84 { 1 } else { 0 })
let response = msi  # perfect separation!

# 1. Try fitting glm("response ~ msi", data, "binomial")
# 2. What happens to the coefficient and its standard error?
# 3. Why is this a problem? (Hint: the MLE doesn't exist)
# 4. How would you handle this in practice?

Key Takeaways

• Logistic regression models binary outcomes by predicting probabilities via the sigmoid function — use it instead of linear regression when Y is 0/1
• Coefficients are on the log-odds scale; exponentiate to get odds ratios (OR): e^β
• An OR of 2.0 means the odds double per unit increase — but this is NOT the same as doubling the probability
• ROC curves show the sensitivity/specificity trade-off across all thresholds; AUC summarizes overall discrimination
• The optimal threshold depends on clinical context: screening favors sensitivity, confirmation favors specificity
• Logistic regression is a special case of the GLM framework — the same approach extends to count data (Poisson), survival data, and more
• Always report odds ratios with confidence intervals, not just p-values
• Watch for separation (a predictor perfectly predicts the outcome), which causes coefficient inflation

What’s Next

Sometimes the outcome isn’t just binary — it’s time-to-event: how long until a patient relapses, how long a cell line survives treatment. Day 17 introduces survival analysis, where censoring makes standard methods fail and Kaplan-Meier curves become essential.

Day 17: Survival Analysis — Time-to-Event Data

The Problem

Dr. Elena Volkov is a cancer genomicist analyzing overall survival in 250 lung adenocarcinoma patients. She has tumor sequencing data and wants to answer: Do TP53-mutant patients survive longer than TP53 wild-type patients?

Her first attempt: compute the mean survival time for each group and run a t-test. But she immediately hits a wall. 40% of patients are still alive at the end of the study. Their survival time is at least 36 months — but she doesn’t know their actual survival time. Dropping these patients biases the analysis (the longest survivors are removed). Using 36 months as their survival time underestimates it.

This is the problem of censoring, and it requires survival analysis.

What Is Censoring?

Right-censoring occurs when the event of interest (death, relapse, progression) has not yet happened at the time of observation. The patient is lost to follow-up, the study ends, or they die of an unrelated cause.

Patients B and C are censored — we know a lower bound on their survival, but not the actual value.

Key insight: Censored observations are NOT missing data. They contain real information (“this patient survived at least X months”). Throwing them away wastes data and biases results. Including them as if the event occurred underestimates survival.

Why Standard Methods Fail

The Kaplan-Meier Estimator

The Kaplan-Meier (KM) estimator is the workhorse of survival analysis. It estimates the survival function S(t) = P(survival > t) as a step function that drops at each observed event.

How it works: At each event time tⱼ:

$$\hat{S}(t) = \prod_{t_j \leq t} \left(1 - \frac{d_j}{n_j}\right)$$

Where:

• dⱼ = number of events at time tⱼ
• nⱼ = number at risk just before tⱼ (alive and not yet censored)

Reading a KM curve:

• Y-axis: proportion surviving (starts at 1.0 = 100%)
• X-axis: time
• Steps down: events (deaths)
• Tick marks: censored observations (patients lost but still contribute until that point)
• Steeper drops: periods of high event rate
• Flat plateaus: stable periods

Median Survival

The median survival is the time at which the KM curve crosses 0.50 — the point where half the patients have experienced the event. It is more robust than the mean because it is less affected by a few extreme values.

Clinical relevance: “Median overall survival was 18 months in the treatment arm versus 12 months in the control arm” is the standard way clinical trials report survival data. It’s the single most important number in oncology clinical trials.

The Log-Rank Test

The log-rank test compares survival curves between groups. It asks: “Is the survival experience significantly different between these groups?”

• H₀: The survival curves are identical
• H₁: The survival curves differ

The test compares observed events to expected events (under H₀) at each time point across the entire follow-up period. It gives most weight to later time points.

Common pitfall: The log-rank test can be non-significant even when curves look different, if the difference is early (then converges) or if curves cross. If hazards are not proportional, consider alternatives like the Wilcoxon (Breslow) test, which gives more weight to early events.

Cox Proportional Hazards Model

The Cox PH model is the regression analog for survival data. It models the hazard (instantaneous event rate) as:

$$h(t|X) = h_0(t) \cdot \exp(\beta_1 X_1 + \beta_2 X_2 + \cdots)$$

Where h₀(t) is the baseline hazard and the exponential term scales it by covariates.

The Hazard Ratio

The key output is the hazard ratio (HR):

$$HR = e^{\beta}$$

Key insight: HR = 2.0 does NOT mean “dies twice as fast” or “survives half as long.” It means that at any given time point, the hazard (instantaneous risk of the event) is 2x higher. The relationship between HR and median survival depends on the shape of the baseline hazard.

Proportional Hazards Assumption

The Cox model assumes that the ratio of hazards between groups is constant over time. If TP53-mutant patients have HR = 1.5, this ratio should hold at 6 months, 12 months, and 24 months.

Violations:

• Curves that cross (treatment effect reverses over time)
• HR that changes with time (e.g., surgery risk is high early, then protective later)
• Delayed treatment effects (immunotherapy often shows late separation)

Adjusting for Confounders

Like multiple regression, Cox models can include multiple predictors:

$$h(t) = h_0(t) \cdot \exp(\beta_1 \cdot \text{TP53} + \beta_2 \cdot \text{Age} + \beta_3 \cdot \text{Stage})$$

This gives the HR for TP53 mutation adjusted for age and stage — a cleaner estimate of its independent effect.

Clinical relevance: A univariate HR for TP53 might be 1.8 (worse survival). But if TP53-mutant tumors also tend to be higher stage, the adjusted HR might drop to 1.3 after controlling for stage. The adjusted HR is what matters for understanding TP53’s independent prognostic value.

Survival Analysis in BioLang

Kaplan-Meier Curves

set_seed(42)
# Lung adenocarcinoma survival by TP53 status
let n = 250

# Simulate TP53 status (40% mutant)
let tp53_mut = rnorm(n, 0, 1) |> map(|x| if x > 0.25 { 0 } else { 1 })

# Simulate survival times (exponential with TP53 effect)
let base_hazard = 0.03  # monthly hazard rate
let survival_time = []
let status = []  # 1 = dead, 0 = censored

for i in 0..n {
    let hazard = base_hazard * if tp53_mut[i] == 1 { 1.6 } else { 1.0 }
    let true_time = -log(rnorm(1, 0.5, 0.2)[0] |> max(0.01)) / hazard
    let censor_time = rnorm(1, 42, 10)[0] |> max(24)

    if true_time < censor_time {
        survival_time = survival_time + [true_time]
        status = status + [1]  # event observed
    } else {
        survival_time = survival_time + [censor_time]
        status = status + [0]  # censored
    }
}

print("Events: {status |> sum} / {n} ({(status |> sum) / n * 100 |> round(1)}%)")
print("Censored: {n - (status |> sum)}")

# Fit Kaplan-Meier for each group
let wt_times = []
let wt_status = []
let mut_times = []
let mut_status = []
for i in 0..n {
    if tp53_mut[i] == 0 {
        wt_times = wt_times + [survival_time[i]]
        wt_status = wt_status + [status[i]]
    } else {
        mut_times = mut_times + [survival_time[i]]
        mut_status = mut_status + [status[i]]
    }
}

let km_wt = kaplan_meier(wt_times, wt_status)
let km_mut = kaplan_meier(mut_times, mut_status)

print("Median survival (TP53 WT): {km_wt.median |> round(1)} months")
print("Median survival (TP53 mut): {km_mut.median |> round(1)} months")

Kaplan-Meier Plot

# Publication-quality survival curves
# Plot KM data using plot()
let km_table = table({
    "time": km_wt.times ++ km_mut.times,
    "survival": km_wt.survival ++ km_mut.survival,
    "group": km_wt.times |> map(|t| "TP53 WT") ++ km_mut.times |> map(|t| "TP53 Mut")
})
plot(km_table, {type: "line", x: "time", y: "survival", color: "group",
    title: "Overall Survival by TP53 Status",
    x_label: "Time (months)", y_label: "Survival Probability"})

Log-Rank Test

# Compare survival curves statistically
# Compute log-rank test from KM outputs
# The test compares observed vs expected events across groups
let km_all = kaplan_meier(survival_time, status)

# Use Cox PH as a proxy for log-rank (equivalent for single binary predictor)
let cox_lr = cox_ph(survival_time, status, [tp53_mut])
print("=== Log-Rank Test (via Cox) ===")
print("p-value: {cox_lr.p_value |> round(4)}")

if cox_lr.p_value < 0.05 {
    print("Significant difference in survival between TP53 groups")
} else {
    print("No significant difference detected (p > 0.05)")
}

Cox Proportional Hazards Model

set_seed(42)
# Simulate additional covariates
let age = rnorm(n, 65, 10)
let stage = rnorm(n, 2.5, 0.8) |> map(|x| max(1, min(4, round(x))))

# Univariate Cox model
let cox_simple = cox_ph(survival_time, status, [tp53_mut])

print("=== Univariate Cox Model ===")
print("TP53 Mutation HR: {cox_simple.hazard_ratio |> round(2)}")
print("  p-value: {cox_simple.p_value |> round(4)}")

# Multivariable Cox model — adjust for age and stage
let cox_adjusted = cox_ph(survival_time, status, [tp53_mut, age, stage])

print("\n=== Multivariable Cox Model ===")
print("Hazard ratios: {cox_adjusted.hazard_ratios}")
print("p-values: {cox_adjusted.p_values}")

# Compare: does TP53 HR change after adjusting for stage?
print("\nTP53 HR unadjusted: {cox_simple.hazard_ratio |> round(2)}")
print("TP53 HR adjusted:   {cox_adjusted.hazard_ratios[0] |> round(2)}")

Forest Plot of Hazard Ratios

# Visualize all HRs from the multivariable model
let hr_data = table({
    "predictor": ["TP53 Mutation", "Age", "Stage"],
    "hr": cox_adjusted.hazard_ratios,
    "p_value": cox_adjusted.p_values
})
forest_plot(hr_data)

# Left of 1 = protective, Right of 1 = harmful

Survival Curve from Cox Model

# Plot adjusted survival curves from KM estimates
let surv_table = table({
    "time": km_wt.times ++ km_mut.times,
    "survival": km_wt.survival ++ km_mut.survival,
    "group": km_wt.times |> map(|t| "TP53 WT") ++ km_mut.times |> map(|t| "TP53 Mut")
})

plot(surv_table, {type: "line", x: "time", y: "survival", color: "group",
    title: "Adjusted Survival Curves from Cox Model",
    x_label: "Time (months)", y_label: "Survival Probability"})

Checking Proportional Hazards

# Check proportional hazards assumption
# If hazard ratios change over time, PH is violated
# Visual check: compare early vs late HR
let early_times = []
let early_status = []
let early_tp53 = []
let late_times = []
let late_status = []
let late_tp53 = []

let midpoint = 24  # months
for i in 0..n {
    if survival_time[i] <= midpoint {
        early_times = early_times + [survival_time[i]]
        early_status = early_status + [status[i]]
        early_tp53 = early_tp53 + [tp53_mut[i]]
    } else {
        late_times = late_times + [survival_time[i]]
        late_status = late_status + [status[i]]
        late_tp53 = late_tp53 + [tp53_mut[i]]
    }
}

print("=== Proportional Hazards Check ===")
print("If HRs differ substantially, PH assumption may be violated")
# If p < 0.05, consider: stratified Cox model, time-varying coefficients,
# or restricted mean survival time (RMST)

Python:

from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test

# Kaplan-Meier
kmf = KaplanMeierFitter()
kmf.fit(time_wt, event_wt, label='TP53 WT')
kmf.plot()
kmf.fit(time_mut, event_mut, label='TP53 Mut')
kmf.plot()

# Log-rank test
results = logrank_test(time_wt, time_mut, event_wt, event_mut)
print(f"p = {results.p_value:.4f}")

# Cox model
cph = CoxPHFitter()
cph.fit(df, duration_col='time', event_col='status')
cph.print_summary()
cph.plot()  # forest plot

# Check PH assumption
cph.check_assumptions(df, show_plots=True)

R:

library(survival)
library(survminer)

# Kaplan-Meier
km_fit <- survfit(Surv(time, status) ~ tp53, data = df)

# Publication-quality KM plot
ggsurvplot(km_fit,
           pval = TRUE, risk.table = TRUE,
           palette = c("#2196F3", "#F44336"))

# Log-rank test
survdiff(Surv(time, status) ~ tp53, data = df)

# Cox model
cox_model <- coxph(Surv(time, status) ~ tp53 + age + stage, data = df)
summary(cox_model)

# Forest plot
ggforest(cox_model)

# Proportional hazards test
cox.zph(cox_model)

Exercises

Exercise 1: Kaplan-Meier and Median Survival

A clinical trial follows 200 breast cancer patients treated with either chemotherapy or chemotherapy + targeted therapy. Compute KM curves and median survival for both arms.

set_seed(42)
let n = 200
let treatment = rnorm(n, 0, 1) |> map(|x| if x > 0 { 1 } else { 0 })  # 0=chemo, 1=chemo+targeted

# Simulate survival data
# Targeted therapy reduces hazard by 30%
# Generate survival times and censoring

# 1. Compute kaplan_meier() for each treatment arm
# 2. Plot KM curves with plot()
# 3. Report median survival for each arm
# 4. What is the absolute improvement in median survival?

Exercise 2: Log-Rank Test with Multiple Groups

Compare survival across four molecular subtypes (Luminal A, Luminal B, HER2+, Triple-negative). Which pairs differ significantly?

set_seed(42)
let n = 300
# Assign subtypes: 0=LumA, 1=LumB, 2=HER2+, 3=TNBC
let subtype = rnorm(n, 0, 1) |> map(|x|
    if x < -0.39 { 0 }
    else if x < 0.25 { 1 }
    else if x < 0.64 { 2 }
    else { 3 })

# Simulate different hazard rates by subtype
# Luminal A: lowest hazard, Triple-neg: highest

# 1. Compute kaplan_meier() for all 4 subtypes
# 2. Fit cox_ph() with subtype as predictor
# 3. Which subtypes have the best and worst prognosis?

Exercise 3: Multivariable Cox Model

Build a Cox model with TP53 status, age, stage, and smoking status. Determine which factors are independently prognostic after adjustment.

let n = 300

# Simulate covariates and survival
# Fit cox_ph() with all 4 predictors
# 1. Report hazard ratios for each predictor
# 2. Which predictors are significant after adjustment?
# 3. Create a forest_plot()
# 4. Does the TP53 hazard ratio change from univariate to multivariable?

Exercise 4: Checking the PH Assumption

Simulate a scenario where the proportional hazards assumption is violated — an immunotherapy that has no effect in the first 3 months but a strong effect afterward (delayed separation). Show that the PH test flags this.

# Simulate two groups:
# Control: constant hazard
# Treatment: same hazard for months 0-3, then 50% reduced hazard

# 1. Plot KM curves — do they cross or show delayed separation?
# 2. Fit cox_ph() and check if HR changes over time
# 3. Is the PH assumption violated?
# 4. What would you do in practice?

Exercise 5: Complete Survival Analysis Pipeline

Perform a full survival analysis: KM curves, log-rank test, univariate Cox for each predictor, multivariable Cox, forest plot. Write a one-paragraph summary of findings.

# Use 400 patients with 5 clinical/genomic variables
# Run the complete pipeline from raw data to clinical interpretation

Key Takeaways

• Censored observations contain real information — survival analysis methods handle them correctly while standard methods cannot
• Kaplan-Meier estimates survival probabilities as a step function; median survival is the primary summary measure
• The log-rank test compares survival curves between groups (the “t-test” of survival analysis)
• Cox proportional hazards regression models the effect of multiple covariates on hazard; the hazard ratio is the key output
• HR = 2.0 means twice the instantaneous risk, NOT twice as fast to die or half the survival time
• Always check the proportional hazards assumption — violations (crossing curves, delayed effects) invalidate the standard Cox model
• Forest plots are the standard visualization for hazard ratios from Cox models
• Adjusted HRs (controlling for confounders) are more clinically meaningful than unadjusted ones

What’s Next

We’ve been analyzing data. But before collecting data, there’s a critical question: How many samples do you need? Day 18 covers experimental design and statistical power — the science of planning studies that can actually detect the effects you’re looking for.

Day 18: Experimental Design and Statistical Power

The Problem

Dr. Ana Reyes is a junior PI writing her first R01 grant. She proposes a study comparing gene expression between psoriatic skin and normal skin using RNA-seq, planning 3 samples per group because “that’s what the lab down the hall used.”

The grant comes back with this reviewer comment:

“The proposed sample size of n=3 per group is inadequate. The applicant provides no power analysis to justify this number. With 3 replicates, the study is severely underpowered to detect anything less than a 4-fold change, which is biologically unrealistic for most genes. We recommend at least 8-10 biological replicates per condition based on published power analyses for RNA-seq DE studies.”

Grant rejected. Six months of proposal writing, wasted — because she didn’t plan the sample size.

Statistical power determines whether your study can actually detect the effect you’re looking for. Getting it wrong wastes time, money, animals, and patient samples.

What Is Statistical Power?

Power analysis sits at the intersection of four interconnected quantities:

The fundamental relationship: Given any three, you can calculate the fourth.

Type I and Type II Errors

Key insight: An underpowered study has a high β — it frequently misses real effects. This doesn’t just waste resources; it can lead to the false conclusion that an effect doesn’t exist, discouraging further research on a real phenomenon.

Why 80% Power?

The convention of 80% power means accepting a 20% chance of missing a real effect. Some contexts demand higher:

Effect Size: The Missing Ingredient

Effect size is the hardest quantity to estimate because it requires knowledge about the biology before doing the experiment. Sources:

Common pitfall: Using an inflated effect size from a small pilot study. Small studies overestimate effects (the “winner’s curse”). If your pilot with n=5 shows d=1.5, the true effect is probably smaller. Be conservative.

Cohen’s d for Two-Group Comparisons

$$d = \frac{|\mu_1 - \mu_2|}{\sigma_{pooled}}$$

Power for Common Designs

Two-Group Comparison (t-test)

The most basic design: compare means between two independent groups.

Required sample size per group (approximate, for α=0.05, power=0.80):

Key insight: Detecting a small effect requires nearly 400 samples per group! This is why GWAS studies need thousands of subjects — individual genetic variants typically have very small effects (d ≈ 0.1-0.2).

Paired Design (Paired t-test)

When the same subjects are measured before and after treatment, pairing removes between-subject variability. Power increases dramatically:

RNA-seq Differential Expression

RNA-seq power depends on additional factors:

Rules of thumb for RNA-seq DE:

• Minimum: 3 biological replicates (detects only >4-fold changes)
• Good: 6-8 replicates (detects 2-fold changes)
• Ideal: 12+ replicates (detects 1.5-fold changes)
• Technical replicates have diminishing returns — invest in biological replicates

Common pitfall: Confusing technical replicates (sequencing the same library twice) with biological replicates (independent biological samples). Only biological replicates give you power to generalize. Ten technical replicates of one sample give you n=1, not n=10.

GWAS

Power Curves

Power curves show how power varies with sample size for different effect sizes. They’re the most informative visualization for study planning — you can see the “sweet spot” where adding more samples gives diminishing returns.

The Cost of Underpowered Studies

An underpowered study is not just a failed study — it’s actively harmful:

1. Waste: Money, time, and irreplaceable biological samples consumed for inconclusive results
2. Publication bias: Only “lucky” underpowered studies (that happen to find p < 0.05) get published, inflating reported effect sizes
3. False negatives: Real treatments or biomarkers get abandoned
4. Ethical cost: Patients enrolled in clinical trials with no realistic chance of detecting a benefit

Clinical relevance: The FDA and EMA require power analyses for all clinical trial protocols. Journal reviewers increasingly require them for observational studies too. “How many samples do you need?” is the first question of good experimental design.

Experimental Design in BioLang

Basic Power Analysis for t-test

set_seed(42)
# How many samples to detect a 2-fold change in gene expression?

# Parameters
let alpha = 0.05
let power_target = 0.80
let effect_size = 0.8   # Cohen's d for ~2-fold change

# Simulate power at different sample sizes
let sample_sizes = [5, 10, 15, 20, 30, 50, 75, 100]
let n_simulations = 1000

print("=== Power Analysis: Two-Sample t-test ===")
print("Effect size (Cohen's d): {effect_size}")
print("Alpha: {alpha}")
print("")
print("n per group    Estimated Power")

for n in sample_sizes {
    let significant = 0

    for sim in 0..n_simulations {
        # Simulate two groups with known effect
        let group1 = rnorm(n, 0, 1)
        let group2 = rnorm(n, effect_size, 1)

        let result = ttest(group1, group2)
        if result.p_value < alpha {
            significant = significant + 1
        }
    }

    let power = significant / n_simulations
    let marker = if power >= power_target { " <-- sufficient" } else { "" }
    print("{n}            {power |> round(3)}{marker}")
}

Power Curves for Different Effect Sizes

set_seed(42)
# Visualize power as a function of sample size
let sample_sizes = [5, 10, 15, 20, 25, 30, 40, 50, 75, 100]
let effect_sizes = [0.3, 0.5, 0.8, 1.2]
let n_sims = 500

let power_curves = {}

for d in effect_sizes {
    let powers = []

    for n in sample_sizes {
        let sig_count = 0
        for s in 0..n_sims {
            let g1 = rnorm(n, 0, 1)
            let g2 = rnorm(n, d, 1)
            if ttest(g1, g2).p_value < 0.05 {
                sig_count = sig_count + 1
            }
        }
        powers = powers + [sig_count / n_sims]
    }

    power_curves["{d}"] = powers
}

# Plot power curves
let curve_table = table({
    "n": sample_sizes,
    "d_0.3": power_curves["0.3"],
    "d_0.5": power_curves["0.5"],
    "d_0.8": power_curves["0.8"],
    "d_1.2": power_curves["1.2"]
})
plot(curve_table, {type: "line", x: "n",
    title: "Power Curves: Two-Sample t-test",
    x_label: "Sample Size per Group", y_label: "Statistical Power"})

RNA-seq Experiment Design

set_seed(42)
# How many biological replicates for RNA-seq DE?

# Simulate RNA-seq-like data
let fold_changes = [1.5, 2.0, 3.0, 4.0]
let replicates = [3, 5, 8, 12, 20]
let n_sims = 200

print("=== RNA-seq Power by Fold Change and Replicates ===")
print("FC       n=3     n=5     n=8    n=12    n=20")

for fc in fold_changes {
    let powers = []

    for n in replicates {
        let detected = 0

        for sim in 0..n_sims {
            # Simulate one gene with known fold change
            let control = rnorm(n, 10, 2)
            let treatment = rnorm(n, 10 * fc, 2 * fc)

            # Log-transform (as in real RNA-seq analysis)
            let log_ctrl = control |> map(|x| log2(max(x, 0.1)))
            let log_treat = treatment |> map(|x| log2(max(x, 0.1)))

            let p = ttest(log_ctrl, log_treat).p_value
            if p < 0.05 { detected = detected + 1 }
        }

        powers = powers + [detected / n_sims]
    }

    print("{fc}     " ++ powers |> map(|p| "{(p * 100) |> round(0)}%") |> join("   "))
}

# Key takeaway: n=3 barely detects 4-fold changes;
# n=8 reliably detects 2-fold changes

Paired vs. Unpaired Design

set_seed(42)
# Show the power advantage of paired designs
let n_sims = 1000
let n = 20
let effect = 0.5  # medium effect
let subject_sd = 2.0  # between-subject variability
let within_sd = 0.5    # within-subject variability

let power_unpaired = 0
let power_paired = 0

for sim in 0..n_sims {
    # Unpaired: independent groups
    let group1 = rnorm(n, 0, subject_sd)
    let group2 = rnorm(n, effect, subject_sd)
    if ttest(group1, group2).p_value < 0.05 {
        power_unpaired = power_unpaired + 1
    }

    # Paired: same subjects, before and after
    let baseline = rnorm(n, 0, subject_sd)
    let after = baseline + effect + rnorm(n, 0, within_sd)
    let diff = after - baseline
    if ttest_one(diff, 0).p_value < 0.05 {
        power_paired = power_paired + 1
    }
}

print("=== Paired vs Unpaired Design (n={n}, d={effect}) ===")
print("Unpaired power: {(power_unpaired / n_sims * 100) |> round(1)}%")
print("Paired power:   {(power_paired / n_sims * 100) |> round(1)}%")
print("Pairing advantage: {((power_paired - power_unpaired) / n_sims * 100) |> round(1)} percentage points")

Multi-Group Design (ANOVA)

set_seed(42)
# Power for detecting differences among 4 treatment groups
let n_sims = 500
let k = 4  # number of groups
let group_means = [0, 0.3, 0.6, 0.9]  # increasing effect
let sample_sizes = [5, 10, 15, 20, 30]

print("=== ANOVA Power (k={k} groups) ===")
for n in sample_sizes {
    let sig = 0

    for sim in 0..n_sims {
        let groups = []
        for i in 0..k {
            groups = groups + [rnorm(n, group_means[i], 1)]
        }

        let result = anova(groups)
        if result.p_value < 0.05 { sig = sig + 1 }
    }

    print("n = {n} per group: power = {(sig / n_sims * 100) |> round(1)}%")
}

Sample Size Recommendation Report

set_seed(42)
# Generate a complete sample size recommendation
let scenarios = [
    { name: "Conservative (d=0.5)", effect: 0.5 },
    { name: "Expected (d=0.8)", effect: 0.8 },
    { name: "Optimistic (d=1.2)", effect: 1.2 }
]

print("=== SAMPLE SIZE RECOMMENDATION REPORT ===")
print("Two-group comparison, alpha=0.05, power=80%")
print("")

for s in scenarios {
    # Find minimum n for 80% power via simulation
    let required_n = 0
    for n in 5..200 {
        let power = 0
        for sim in 0..500 {
            let g1 = rnorm(n, 0, 1)
            let g2 = rnorm(n, s.effect, 1)
            if ttest(g1, g2).p_value < 0.05 { power = power + 1 }
        }
        if power / 500 >= 0.80 {
            required_n = n
            break
        }
    }

    print("{s.name}: n = {required_n}/group")
}

print("")
print("Recommendation: Plan for the CONSERVATIVE")
print("estimate + 10-20% for dropout/QC failures")

Python:

from scipy.stats import norm
from statsmodels.stats.power import TTestIndPower, TTestPower

# Power analysis for two-sample t-test
analysis = TTestIndPower()

# Required sample size
n = analysis.solve_power(effect_size=0.8, alpha=0.05, power=0.80)
print(f"Required n per group: {n:.0f}")

# Power curve
import matplotlib.pyplot as plt
fig = analysis.plot_power(
    dep_var='nobs', nobs=range(5, 100),
    effect_size=[0.3, 0.5, 0.8, 1.2])

# Simulation-based power
import numpy as np
from scipy.stats import ttest_ind

def simulate_power(n, d, n_sim=1000):
    sig = sum(ttest_ind(np.random.normal(0, 1, n),
                        np.random.normal(d, 1, n)).pvalue < 0.05
              for _ in range(n_sim))
    return sig / n_sim

R:

# Power analysis for two-sample t-test
power.t.test(d = 0.8, sig.level = 0.05, power = 0.80)

# Power curve
library(pwr)
pwr.t.test(d = 0.8, sig.level = 0.05, power = 0.80)

# RNA-seq specific power
library(RNASeqPower)
rnapower(depth = 20e6, cv = 0.4, effect = 2,
         alpha = 0.05, power = 0.8)

# Simulation-based
library(simr)
# simr provides power simulation for mixed models

Exercises

Exercise 1: Power for Your Study

You’re planning a study comparing tumor mutation burden between immunotherapy responders and non-responders. Pilot data suggests d ≈ 0.6 with SD = 5 mutations/Mb. How many patients per group do you need for 80% power?

# 1. Simulate power at n = 10, 20, 30, 50, 75, 100
# 2. Find the minimum n for 80% power
# 3. Add 15% for anticipated dropout
# 4. Create a power curve plot