I was looking at some analysis from a colleague recently where they tested 50 different variables to see if they correlated with customer churn. They came back thrilled, reporting that seven of those variables were "statistically significant" (p < 0.05). My immediate reaction wasn't excitement—it was skepticism. When you run 50 tests, the laws of probability dictate that you'll find some "winners" even if the data is complete garbage.

The "Wait, that can't be right" moment

Let's simulate this. I'll create a scenario where we are testing 20 different "wellness supplements" to see if they lower blood pressure. In reality, I'm going to make it so that 19 of these supplements do absolutely nothing, and only one actually works. We should only find one significant result, right?

# Setup: 20 supplements, 10 people per group
set.seed(42)
results <- replicate(20, {
  # Group A: Placebo, Group B: Supplement
  # Most are just random noise (mean=120 for both)
  group_a <- rnorm(10, 120, 5)
  group_b <- rnorm(10, 120, 5)
  t.test(group_a, group_b)$p.value
})

# Let's force the 20th one to actually be significant
results[20] <- 0.001 

# How many look "significant" at the 0.05 level?
sum(results < 0.05)

When I run this, I often see 2 or 3 "significant" results. I know for a fact that only one supplement actually works, yet R is telling me others are significant too. This is the "Multiple Comparisons Problem." If you roll a die 20 times, you're probably going to hit a six eventually; that doesn't mean the die is magic, it just means you tried a lot of times.

Applying the Sledgehammer: Bonferroni

The most intuitive way to fix this is to be way stricter. If I'm doing 20 tests, I'll just divide my alpha (0.05) by 20. This is the Bonferroni correction. It's the "sledgehammer" approach—it kills false positives, but it might kill your actual discoveries too.

# Using the built-in p.adjust function
bonf_p <- p.adjust(results, method = "bonferroni")
sum(bonf_p < 0.05)

Now we're likely back to just 1 significant result. It worked! But here's the catch: Bonferroni is incredibly conservative. If I had 1,000 variables (like in genomic data), my adjusted p-value would have to be 0.00005 to be considered significant. I'd end up throwing away real, valuable signals because the bar was set impossibly high.

Finding a Middle Ground with FDR

I usually prefer the Benjamini-Hochberg (BH) method. Instead of trying to ensure there are zero false positives (which is what Bonferroni does), BH controls the False Discovery Rate. It basically says, "I'm okay if 5% of my 'significant' results are fake, as long as the rest are real."

Let's see how the BH method handles our supplement data compared to the others:

# Comparing methods
comparison <- data.frame(
  Raw = results,
  Bonferroni = p.adjust(results, method = "bonferroni"),
  BH = p.adjust(results, method = "BH")
)

head(comparison)

If you look at the BH column, the p-values are higher than the raw ones, but lower than the Bonferroni ones. It's a sliding scale. In a professional pipeline, this is almost always the right choice when you're doing exploratory data analysis. It gives you the power to find effects without being fooled by the sheer volume of tests you're running.

Which one do you actually use?

I've found a simple rule of thumb over the years: If the cost of a False Positive is catastrophic (e.g., approving a drug that might be toxic), use Bonferroni. If you're in the "discovery" phase—trying to find interesting leads to investigate further—use BH. Just remember that any time you see a p-value in a paper or a report, your first question should be: "How many other tests did they run before they found this one?"

📋 Practical Task

Filtering Differentially Expressed Genes

You are analyzing a transcriptomics dataset. You have a vector of p-values representing 100 different genes. You suspect that only a few are actually differentially expressed, but the rest are noise.

Your task:

1. Create a vector of 100 p-values. Make 95 of them random values between 0.01 and 1.0, and 5 of them very small (e.g., between 0.0001 and 0.001).
2. Apply the Benjamini-Hochberg (BH) correction using p.adjust().
3. Create a logical vector that identifies which genes remain significant at an adjusted alpha level of 0.05.
4. Print the number of genes that passed the raw 0.05 threshold versus the number that passed the BH-adjusted 0.05 threshold.