2
$\begingroup$

I referred to this question, but found no solution to my problem.

I want to simulate a simple case of multiple comparison problem, where family-wise alpha error increases.

Question 1

First example includes only a single test:

library (tidyverse)

# 100 datasets with 100 participants, means and sd
set.seed(1)
A <- rerun(.n = 1000, rnorm(n = 100, mean = 0, sd = 1))
B <- rerun(.n = 1000, rnorm(n = 100, mean = 0, sd = 1))
C <- rerun(.n = 1000, rnorm(n = 100, mean = 0, sd = 1))

test1 <- map2_dbl(A, B, ~t.test(.x, .y)$p.value)
sum(test1<=0.05 )/length(test1) # 0.044

However, when i try to simulate a case where I compare A-B, A-C and B-C with t-tests, I don't get the expected amount of significant p-values $\alpha_F=1-(1-0.05)^c$, where c is the number of comparisons (here - 3). Family-wise alpha should be $\alpha_F=0.1426$.

Here's the code:

triple_t <- function(x,y,z) {
  t1 <- t.test(x,y)
  t2 <- t.test(x,z)
  t3 <- t.test(y,z)
  return(c(t1$p.value, t2$p.value, t3$p.value))
}

test2 <- pmap(list(A, B, C), triple_t) %>% flatten_dbl()
sum(test2<=0.05 )/length(test2) # 0.047

What am I missing here?

Question 2

If we have a situation with 6 different datasets, all drawn from the same population

set.seed(1)
A <- rerun(.n = 1000, rnorm(n = 100, mean = 0, sd = 1))
B <- rerun(.n = 1000, rnorm(n = 100, mean = 0, sd = 1))
C <- rerun(.n = 1000, rnorm(n = 100, mean = 0, sd = 1))
D <- rerun(.n = 1000, rnorm(n = 100, mean = 0, sd = 1))
E <- rerun(.n = 1000, rnorm(n = 100, mean = 0, sd = 1))
F <- rerun(.n = 1000, rnorm(n = 100, mean = 0, sd = 1))

would you expect the same family-wise error probability to increase if you conduct comparisons only on independent samples (A-B, C-D, E-F)?

triple_talt <- function(a,b,c,d,e,f,g) {
  t1 <- t.test(a,b)
  t2 <- t.test(c,d)
  t3 <- t.test(e,f)
  return(c(t1$p.value, t2$p.value, t3$p.value))
}

test4 <- pmap(list(A, B, C, D, E, F), triple_talt) 
$\endgroup$
0

1 Answer 1

2
$\begingroup$

The problem seems to be with how you are interpreting the family-wise error rate. There are three tests being conducted simultaneously. Let $E_i$ be the event that hypothesis $i$ is rejected, $i=1,2,3$. In your simulations, all null hypotheses are true, so that $$P(E_i) = 0.05, \ i=1,2,3.$$

What your code is doing.

You are running $M=1000$ simulations, counting the total number of rejections and dividing by $3M$. You are essentially obtaining a Monte Carlo estimate of the quantity $$\frac{P(E_1)+P(E_2)+P(E_3)}{3},$$ which is equal to $0.05$ by construction.

What you should be doing.

The family-wise type I error rate refers to the quantity $$P(E_1 \cup E_2 \cup E_3).$$ Since all hypotheses are true, this is the probability that there is at least one false positive, each time data is collected. The following adjustment to your code will accomplish this.

test2 <- pmap(list(A, B, C), triple_t) %>% lapply(function(z) min(z) < 0.05)
mean(unlist(test2))

Running this adjusted code gives a family-wise type I error of $0.111$.

The theoretical error rate

You might be wondering why this error rate, $0.111$, differs from the theoretical error rate $0.1426$. The formula you cite can be derived as \begin{align*} P(E_1 \cup E_2 \cup E_3) &= 1 - P(E_1^c \cap E_2^c \cap E^c) \\ &= 1 - P(E_1^c)^3 \\ &= 1 - (1-0.05)^3 \end{align*}

This error bound, known as a Sidak correction, is a valid upper bound on the true error rate when the tests are positively dependent. The formula will only be exact if the tests are independent. The tests here are not independent due to the shared variables between tests. For instance, assume that we lack evidence to conclude that $\mu_A \neq \mu_B$ and we also lack evidence to conclude that $\mu_A \neq \mu_C$. Conditional on this information, the probability that we reach the conclusion $\mu_B \neq \mu_C$ is now less than $0.05$. If each of the three tests was based on independent variables, i.e. $A_1-A_2$, $B_1-B_2$, $C_1-C_2$, you would achieve the error rate $0.1426$ exactly.

$\endgroup$
2
  • $\begingroup$ Thank you for the detailed answer - could you please take a look at Question 2 that is related to the same topic? $\endgroup$
    – User33268
    Feb 25, 2020 at 13:47
  • $\begingroup$ @User33268 Much of that question is answered in the final paragraph of my answer. To clarify: In question 1, the tests are "positively dependent" so the Sidak FWER gives an upper bound. So in question 2, we would expect the error rate to increase (to exactly $0.1426$ if you do enough simulations) because the Sidak FWER is exact for independence tests. If the tests were negatively dependent (less common in practice, but still possible), then the error rate would decrease when you make the tests independent. $\endgroup$
    – knrumsey
    Feb 25, 2020 at 18:41

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.