simulating a multiple comparisons problem using R and bonferroni correction

Question

I am actually working on an exercise for my students dealing with multiple comparisons. However, I was having trouble trying to simulate some data in R to demonstrate the problem with multiple comparisons. I was hoping someone might have some R code or at least give some guidance on how to generate a dataset in R to demonstrate the multiple comparisons problem. So first I would like to generate some data (rnorm(), runif, etc.). Then I would like to show how if I pairwise compare the parameters for 2 levels of a multiple level predictor, I am more likely to get a significant result than the standard alpha level of 0.5, or 95% confidence.

Thanks.

Answer 1

Let us simulate $100$ realisations of $N(0, 1)$ and test the null hypothesis $H_0: \mu=0$ with the t-test. If this is done many times, $H_0$ should be rejected approximately $5\%$ of the times.

pval <- rep(NA, 1e4)
for(k in 1:1e4)
{
  x <- rnorm(100)
  pval[k] <- t.test(x=x, mu=0)$p.value
}

> mean(pval < 0.05)
[1] 0.0495

$$ $$

If we now simulate $5$ random variables and test the null hypothesis that all means are simultaneously $0$, then the probability of at least one significant result is larger than $5\%$; actually it is $1 - (1 - 0.05)^5 = 0.226$.

n <- 5
pval <- matrix(NA, nrow=1e4, ncol=n)
for(k in 1:1e4)
{
  X <- replicate(n=n, expr=rnorm(100))
  pval[k, ] <- apply(X=X, MARGIN=2, FUN=function(x){
    t.test(x=x, mu=0)$p.value
  })  
}

> mean(sapply(X=1:nrow(pval), FUN=function(i){
+   any(pval[i, ] < 0.05)
+ }))
[1] 0.2203

$$ $$

But if you use the Bonferroni correction, then

> mean(sapply(X=1:nrow(pval), FUN=function(i){
+   any(pval[i, ] < 0.05 / n)
+ }))
[1] 0.0487

Answer 2

You are in Luck I had the same question yesterday myself. I build a function that lets you visualize the issue.

# Say you have a range of tests and you would want to see the probabilty of multiplicity
ntests <- (1:50)

  multicomp.issue <- function(a,k) { # a is my alpha level, k is how many different tests i am doing
    b <- rep(NA,length(k))
    for (i in k) { 
      b[i] <-1-(1-a)^k[i] # multicomp issue formula
    } 

    plot(b, ylim = 0:1, 
         ylab = "Probability",
         main ="Multiplicity error", type = "p", pch = 15)
    abline(h=0.5, col = "red") # 50% 
    abline(h=0.05, col = "green") # what people usually use as significance level 

    print("The family wise error ratre is equal to") # this shows you the numbers in the console. 
    print(b)}
multicomp.issue(0.05, tries)  # Test

When you run that code in your R you should get the following plot.
Of course I advise you play with different a values and different amounts of tests.
I hope this is usefull to you.