3
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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
| cite | improve this answer | |
$\endgroup$
1
$\begingroup$

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.

enter image description here

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.