Multiple hypothesis testing

Question

Suppose I have time series data for $(X,Y_i)$, independent and dependent variables, respectively, which vary over time $t=1,...,1000$. I have data on $Y_i$ for 100 individuals/cities/etc $i=1,...100$. So, each of these 100 $Y$'s is a time series and has 1000 data points. $X$ is common across all individuals and also varies over time.

I run 100 regressions of $Y_i$ on $X$, one for each $i$. I find that 1 one of those regressions is significant at, say, $\alpha=10\%$ level.

Since the probability under a binomial with m=100 trials and probability of success ($\alpha = 10\%$) of having at least one success (one false positive) is almost 1, this does not mean much. The probability that there is at least one false positive in one of my 100 regressions is almost 1.

A common way to go about this is to use a Bonferroni kind of correction that would demand the use of a $\alpha / m$ significance level.

However, the probability of finding more than 15 false positives is $7.25\%$, which is actually lower than $10\%$. Therefore, if I would have found more than 15 significant coefficients across my 100 regressions, I could be relatively confident that I was not facing a chance result.

Is this reasoning justified or not? In fact, in the sense of decreasing the type I error, my reasoning would suggest that I would be more confident that I am not getting false positives if I found 15/100 significant coefficients at the $10\%$ level compared to the situation in which I get 1/1 signifcant coefficient in one single regression at the $10\%$ level. In such a case, I would be inclined to say that a Bonferroni correction would not make sense (it would be unncessary).

A bit more detail as requested:

I run 100 time series regressions of the form (one for each $i$):

$Y_{i,t} = \alpha + \beta_iX_t$, $i=1,2,3,...,100$

$Y_{i,t}$ is happiness of country $i$ in month $t$. $X_{t}$ is world GDP.

Note that this is not a panel stacking countries. This is one time series regression per country. For instance, for China:

$Y_{china,t} = \alpha + \beta_{china}X_t$

I reject the null that $\beta_i= 0$ at the 10% level for some countries $i$ and not others.

The probability that in 100 trials I reject at least one by chance, is 99% under a binominal distribution with probability of success (false positive) of $10\%$ and 100 trials. Therefore, if I would find 1 out of 100 significant coefficients, the probability of type I error is almost 1: out of 100 tries, I'd expect almost surely at least one false positive anyways. So I should not reject the null that $\beta_i= 0$ for that country, even though it passes the individual significance test at the 10\% level.

Now, suppose that instead of 1 rejection, I am actually able to reject the null for 15 countries. The probability that in 100 trials I reject at least 15 by chance, is $7.25%$ under a binominal distribution with probability of success (false positive) of 10% and 100 trials. Therefore, I am tempted to conclude that, given that I am willing to live with a type I error chance of 10% (just assume this), I can in fact reject the null $\beta_i= 0$ for those 15 countries and indeed say that I have evidence that x is associated with y.

Therefore, what I am saying, is that I do not understand why I would ever apply a Bonferroni or FDR correction in this setting. If i I have only one significant coefficient out of 100, I cannot trust it because I'd expect it with probability 1. I would need some kind of correction (like Bonferroni)! But if I have 15 significant coefficients, I can actually trust it more since the probability of observing that by chance is less than 10\%. I do not need need any correction. Take the extreme case: instead of 15, all 100 coefficients are individually significant at the 10\% level. The prob that at least 100 coefficients out of 100 are significant (under binomial) is almost 0%. I couldn't possibly be more confident that my results are not driven by type I error. Why would I want to adjust or correct for anything?

My question is simply: is this reasoning wrong or correct? I claim I do not need a correction. Yes, or no? If no, where is the mistake in my reasoning. Thanks!

New add-on as reply to greg (no space in the comment):

Thank you Greg! Awesome insightful explanation. I am just still confused on why it is such a big problem to assume that all null hypotheses are true. Isn't that what we do when we compute p-values: the probablity of observing a value more extreme than the one we estimated, given that the null is true? My idea is to basically derive a joint p-value by defining a probability of observing 15 values more extreme than the 15 ones I estimated across 100 regressions, given that the 100 null hypotheses are true (the worst case scenario in which everything is insignificant).

If I understand your explanation correctly, the problem with my approach would be that, focusing on type I error only, for instance, maybe 7 null hypotheses are true and 8 are false. I reject 15, so I have 7 false positives. The probablity of having at least 7 false positives out of 92 coefficients, is 82\%, which is much larger than the desired 10\% level I would like. Is this correct?

Answer 1

You are thinking along the correct lines, but it looks like you have more to learn as well. This is a topic that you can spend a lot of time on (much more that can really be answered in this type of forum). I would suggest the Wikipedia page https://en.wikipedia.org/wiki/Multiple_comparisons_problem as a place to start, then follow the links and read the references from there.

The Bonferroni is known to be conservative in most cases (it over adjusts so that the real probability of getting any type I errors is less than the target rate as you computed with the 7.25% compared to 10%). The False Discovery Rate and others are becoming more popular, but require thinking about it a little differently.

Most of the standard corrections just increase the p-value or make the confidence intervals wider to take into account the extra variability from the multiple tests. But these do not do anything to take into account that the most "significant" estimates are likely to be biased away from the null by chance.

There are penalized regression methods (ridge, lasso, etc.) that bias the extreme estimates towards 0 (the null) in an attempt to counter the multiple comparison bias, but these do not provide nice simple p-values or confidence intervals (but can be very good for prediction modelling and general exploration).

Bayesian hierarchical models provide another option for inference where you fit all the models together and let them learn from each other. This also pulls the extreme estimates towards the center for a more meaningful model.

Edit for expanded question.

The way that your expanded question is worded indicates that you are still mis-understanding some fundamentals of frequentist hypothesis testing.

All your calculations are conditional on the number of hypotheses that are actually true or false (just unobserved). You talk about the probability of rejecting at least one hypothesis out of 100 by chance as being 99%, but that is only if all the null hypotheses are true. If they are all false then the probability of rejecting any by chance alone is 0%. If half are true and half are false then the probability of rejecting by chance is between 0% and 99%.

You are partially correct in that if you reject one null out of 100 that that can easily be explained by that (and others) being true and just chance gave you significance. But it is also possible that that null is false and you correctly rejected it. It is even possible that many of the other nulls are also false and by not rejecting them you made type II errors. So basically in this case there are many possibilities and you have not ruled any of them out.

You talk about the probability of having made a type I error conditional on rejecting hypotheses, but you can't do that under frequentist theory. Each null hypothesis is either true or false (just unknown which) so the probability of having made a type I error is either 0% or 100% (just unknown which). If 15 of the 100 were significant this could be because all 15 nulls were false and you made correct decisions, all 15 nulls were true and you made a type I error (this is unlikely due to chance, but still possible), or some subset of the 15 are true and the rest false. Part of the idea of the False Discovery Rate is to say I know that I will falsely reject some nulls (and possibly correctly reject some in the process as well) I just want to limit how many false positives I will find (on average). If you have the false discovery rate set at 10% and reject 15 out of 100 nulls, this would imply that 10 of those rejections are false positives and 5 are true positives (but this is on average, so it could be 9 and 6, or 11 and 4, or any other combination). The usual follow-up question is then "which 5 are real?" and the answer is "You don't know".

Now if your goal is just to say "at least 1 slope is not equal to 0" and you don't care specifically which one, then your binomial calculations can work and you can say 15 out of 100 is unlikely if all the nulls were true, so you reject the composite null of all 100 being true. This is a form of multiple comparison adjustment (based on different assumptions than the Bonferroni, but still legitimate if you are happy with the assumptions, this is one form of meta-analysis). Your final question is that you think you don't need a correction, but you made a correction, just different from the Bonferroni (you don't need the Bonferroni in addition to what you did).

But note that your slope estimates (conditional on rejecting the null) are likely biased away from 0. You can see this for yourself with some simulations.

Here is some simple R code:

set.seed(1)
x <- 1:25
y <- replicate(100, y <- 2 + 1*x + rnorm(length(x),0,65), simplify=FALSE)

betas <- sapply(y, function(yy) {
  coef(summary(lm(yy ~ x)))[2,c(1,4)]
})

# how many significant at alpha=0.1
sum(betas[2,] <= 0.1)

# all betas (true beta = 1)
hist(betas[1,])
summary(betas[1,])

# Just the "significant" ones (remember true beta=1)
hist(betas[1, betas[2,] <= 0.1])
summary(betas[1, betas[2,] <= 0.1])

In this case all the slopes are 1, 100 regressions are fit and 16 have p-values less than 0.1. If you look at the summary of all the slopes they are centered close to 1, but if you only look at the "significant" ones, then they are all greater than 2, quite a bit of bias. You can change the seed and other pieces of this simulation and try it multiple times for yourself. The bias is not a problem if you only want to say at least one slope is not 0, but as soon as you start to try to predict new values or interpret any of the "significant" slopes, this will be a problem.

You talk about the probability of the null hypothesis being false (or true), then you are starting to think more along Bayesian than Frequentist lines. In that case it is better to just do a full Bayesian approach (done properly this can also help with the bias in the slope estimates). But you should do the Bayesian approach properly with a meaningful prior (this takes some thought).

Edit 2

Maybe these 2 references will help:

https://www.sciencedirect.com/science/article/pii/S0140673605664616 https://www.sciencedirect.com/science/article/pii/S0140673605665166