I have 10 models analysing pupal developmental rates (DR) and developmental times (DT) depending on a number of factors:

m1 <- lm(DR ~ population+sex+temperature+weight+population*temperature)
m2 <- lm(DR ~ population+sex+population*sex)
m3 <- lm(DR ~ population+sex+temperature+population*temperature)
m4 <- gls(DR ~ population+sex, ...)
m5 <- lm(DR ~ population+sex+temperature+weight+population*temperature)
m6 <- glm(DT ~ sex+temperature+weight, gamma)

... and so on

The models refer to datasets of two different species, two different years and different phases of pupal development.

One of the reviewers asked for control of the experiment-wide false discovery rate, which is not so common in the field of ecology. In many of the examples from genetics and medical studies one testing results in only one p-value, which is then adjusted using the FDR procedures. How does it work in my case? Extracting p-values of only, say, the population variable and use them? And then the same separately for sex and so on? Or use all the p-values from each of the models as one vector?

If so, what shall I do in cases when the factor has already been excluded from the final model – skip them and use only the ones left? Add them back to the model?


2 Answers 2


@Dian breathe easy, it's pretty much not too difficult. So let's work from familiar territory to false discovery rate (FDR).

First, I see that you have a bunch of outcomes, with a varying number of predictors. Someone who is more familiar with multivariate regression (i.e. multiple dependent variables, assuming possible correlations between errors of different models) will have to speak to whether your modeling approach is the best one. Let's take it as given.

Each of your models will produce some number of $p$-values (incidentally I am an epidemiologist, and have absolutely no idea what you mean about "only one $p$-value." If that were true, it would change the nature of my work and that of my colleagues considerably. :). You could go ahead and test your hypotheses about individual effects separately using these $p$-values.

Unfortunately, hypothesis testing is like the lottery (the more you play, the more your chance to "win"), so if you want to go into each hypothesis test assuming that the null hypothesis is true, then you are in trouble, because $\alpha$ (your willingness to make/probability of making a false rejection of a true null hypothesis) only applies to a single test.

You may have heard of "the Bonferroni correction/adjustment", where you try to solve this conundrum by multiplying your $p$-values by the total number of null hypotheses you are testing (let's call that number of tests $m$). You are effectively trying to redefine $\alpha$ as a family-wise error rate (FWER), or the probability of making at least one false rejection out of a family of tests, assuming all null hypotheses are true. Alternatively, and equivalently, you can think about the Bonferroni adjustment as dividing $\alpha$ by $m$ (or $\alpha/2$ by $m$ if you are performing two-tailed tests, which in all likelihood you are in a regression context). We get these two alternatives because basing a rejection decision on $p \le \frac{\alpha/2}{m}$ is equivalent to $mp \le \frac{\alpha}{2}$.

Of course, the Bonferroni technique is a blunt hammer. It positively hemorrhages statistical power. $\overset{_{\vee}}{\mathrm{S}}\mathrm{idák}$ got a smidge more statistical power, by altering the adjustment of the $p$-value to $1-(1-p)^{m}$. Holm improved upon both Bonferroni and $\overset{_{\vee}}{\mathrm{S}}\mathrm{idák}$ adjustments by creating a stepwise adjustment. The step-up procedure for the Bonferroni adjustemnt:

  1. Compute the exact $p$-value for each test.
  2. Order the $p$-values from smallest to largest.
  3. For the first test, adjust the $p$-value to be $pm$; and generally:

    For the i$^{\text{th}}$ test, adjust the $p$-value to be $p(m–(i–1))$.

  4. Using Holm’s method, for all tests following the first test for which we fail to reject H$_{0}$ we will also fail to reject the null hypothesis.

The Holm-$\overset{_{\vee}}{\mathrm{S}}\mathrm{idák}$ adjustment is similar, but you would adjust each $p$-value using $1-(1-p)^{m-(i-1)}$.

Some folks, most notably Benjamini and Hochberg (1995), were not comfortable with the world view implied by the assumption that all null hypotheses are true within a stepwise procedure. Surely, they reasoned, if you make an adjustment and reject a single hypothesis, that must imply that a better assumption would be that the remaining $m-1$ hypotheses have a lower probability of all null hypotheses being true? Also, science in general does not assume that there are no relationships in the world: quite the opposite, in fact. Enter the FDR which progressively assumes that rejection probabilities must increase if previous hypotheses were rejected after adjustment. Here's the step-down procedure they proposed:

  1. Compute the exact $p$-value for each test.
  2. Order the $p$-values from largest to smallest (step-down!).
  3. For the first test ($i=1$), adjust the $p$-value to be $\frac{pm}{m-(1-1)} = p$.
  4. For the i$^{\text{th}}$ test, adjust the $p$-value to be $\frac{pm}{m-(i–1)}$.
  5. Using Benjamini & Hochberg’s method, we reject all tests including and following the first test for which we reject the null hypothesis.

We often term $p$-values that have been adjusted this way $q$-values.

The advantages of this FDR adjustment include (1) more statistical power, especially for large $m$, and (2) easy integration of additional tests/$p$-values (say, adding $p$-values from an additional regression model) in a manner which leaves the inferences from the first FDR adjustment unchanged.

Update: All these FWER procedures, and the FDR procedure I just described can produce adjusted $p$-values that are greater than one. When reporting adjusted $p$-values, these are typically reported as $p=1$, $p>.999$, $p=$not reject or something along those lines.


Benjamini, Y. and Hochberg, Y. (1995). Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society. Series B (Methodological), 57(1):289–300.

  • 1
    $\begingroup$ The Bonferroni correction doesn't get the maths wrong: it's a conservative correction for tests that aren't independent. Sidak's is exact, but only for independent tests. $\endgroup$ May 9, 2014 at 16:11
  • $\begingroup$ Fair enough, got confused for a moment. Edited. $\endgroup$
    – Alexis
    May 9, 2014 at 16:30
  • $\begingroup$ @Alexis: thank you so much for your detailed answer. Bonferroni correction is really more commonly used by us, while with the FDR I first encountered only now. I think I'd use BH adjustment, but in the step "Order the $p$-values from largest to smallest" I wonder in my case which $p$-values I should order? In the model output I've got a p-value for population variable, one for sex variable and so on. In the reference paper they consider the example where the treatment groups were compared with chisq and t tests, thus got one test - one pvalue. $\endgroup$
    – Dian
    May 10, 2014 at 7:28
  • 1
    $\begingroup$ @Dian... not at all. Some comments: (1) I don't have a dog in this race, but if a reviewer explicitly asked for FDR, I would go ahead and do FDR :); (2) these methods are implemented in most statistical software, so if you ain't doing by hand, consider that (e.g. p.adjust() in R, qqvalue in Stata, etc.), (3) you have a $p$-value for each variable and interaction term in each model (I counted 20-something in your examples above), you can also add each whole model $p$-value if reporting; (4) if my answer is useful, I hope you will consider clicking the check-mark at the top left of it. $\endgroup$
    – Alexis
    May 10, 2014 at 13:26
  • 1
    $\begingroup$ Clarifying that my (3) in the above comment means you pile together the $p$-values from all the tests you are conducting. $\endgroup$
    – Alexis
    May 10, 2014 at 13:39

What you're doing is a type of step-wise regression, so FDR in this context serves to guide the model selection: the p-value would be that of the goodness-of-fit (coefficient of determination), not of any individual explanatory variable.

  • $\begingroup$ This is a very important question because most analysis are performed with a multivariate model and backwards stepwise selection. And do we need to apply FDR corrections if we are using regularization (lasso)? $\endgroup$
    – skan
    Jul 9, 2023 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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