Are tests of correlation to multilevel factors correlated?

Problem

Suppose we have an obervation and we want to test whether there is a treatment effect from 4 exclusive treatments. The null hypothesis is that the observations are random and independent of the treatment. I will generate 5000 random datasets with the null hypothesis and see whether the tests for each treatment are independent.

pvals = c()
betas = c()
for(i in c(1:4000)) {
  df = data.frame(obs=rnorm(100),group=as.factor(rep(c(0:3),25)))
  aa = summary(lm(formula = obs ~ group, data = df))
  pvals = cbind(pvals,aa$coefficients[,4])
  betas = cbind(betas,aa$coefficients[,1])
}

Correlations

Let us find the correlation of the p-values:

cor(t(pvals))

##             (Intercept)    group1    group2    group3
## (Intercept)   1.0000000 0.3770484 0.3641579 0.3729794
## group1        0.3770484 1.0000000 0.1559180 0.1783368
## group2        0.3641579 0.1559180 1.0000000 0.1666983
## group3        0.3729794 0.1783368 0.1666983 1.0000000

Correlation of the \(\beta\) estimates

cor(t(betas))

##             (Intercept)     group1     group2     group3
## (Intercept)   1.0000000 -0.7028761 -0.7058219 -0.7130741
## group1       -0.7028761  1.0000000  0.4917052  0.5113477
## group2       -0.7058219  0.4917052  1.0000000  0.5140349
## group3       -0.7130741  0.5113477  0.5140349  1.0000000

Implications for testing a composite hypothesis

What is the false discovery rate for testing the effect of 1 treatment using a p-value of 0.05:

mean(pvals[2,]<0.05)

## [1] 0.05

Which is what we would expect. Now what if I am interested in the hypothesis that treatment 1 AND 2 are significant?

mean( (pvals[2,] < 0.05)&(pvals[3,]<0.05) )

## [1] 0.00875

Well, if they were independent, I would expect 0.025 of the time to get p-values less than 0.05 for both cases. But I actually get a number higher than that. So its a little more liberal than independence. But its less than 0.05 so we’ll still be conservative if specified the same threshold.