tl;dr In a random-slopes model, how should one adjust for multiple comparisons when performing inference on the group-specific slopes (the BLUPs)?

Note 1: Bretz et al, the R package 'multcomp', and several other questions on this site deal with multiple comparisons in the context of the fixed effects in mixed-effect models. This question is about the random effects.

Note 2: this question is easiest to ask using the well-developed frequentist vocabulary surrounding surrounding false-discovery rates and multiple comparisons. However, I am equally interested in answers that can provide some Bayesian perspective on this problem: how to temper interpretation of credible intervals in light of the multiple comparisons issue?

Note 3: I've edited the question substantially thanks to helpful comments from Alexis.


Suppose I fit a well-specified random slopes model to data with N groups. I wish to perform inference on whether these slopes (the BLUPs) differ from zero while controlling the false discovery rate. Is it possible to construct a more powerful test than we would achieve using standard p-value adjustments?

Several lines of reasoning suggest that it should be possible, at least in some cases. Below, I present some simulation results suggesting that this is the case, but interpreting the simulations is complicated, for reasons that I'll discuss. First, some conceptual ramblings:

  1. It's useful to distinguish two cases: one where the true value for every group-level slope is zero, and another where it is not. In the first case (every slope is zero) lme4 (and possible other mixed modeling software?) will likely conclude that the random-effect variance is zero, the model collapses to a global-slope model, and the multiple comparisons problem disappears. Thus, counterintuitively, we will face primarily the multiple comparisons problem when there is good evidence in the data for variation in the BLUPs. This, in turn, tends to happen only when some of the true effect sizes are nonzero.

  2. If the random-effects mean is statistically indistinguishable from zero, perhaps one could perform a likelihood ratio test to determine whether the random effect belongs in the model at all. If we can reject the null hypothesis that all groups are equal (i.e. that the random effect is superfluous), then it must be the case that at least some of the group-level slopes are nonzero. And then we might follow-up with some sort of post-hoc test?

  3. The random-slopes model provides a shrinkage estimator for the slopes. If the true random-effect mean (i.e. the fixed effect in lme4 parlance) is zero, this shrinkage should tend to make it harder to reject a true null. Here, Gelman explores this topic in a Bayesian context. If, on the other hand, the true random-effect mean is nonzero, this should make it relatively frequent to reject a true null, because groups with no true effect will tend to get pulled towards the overall mean.

  4. If we want to study the issue using simulation, there's a bit of a problem. We need to inject a large number of true nulls into a single model in order to study the behavior. At the same time, we need to inject some groups with true effects to force the model to estimate a nonzero random effects variance. When we do both of these things, the true distribution of the effect sizes is no longer normal, which is a form of lack-of-fit in the model.

The above caveats (especially in point #4) notwithstanding, I've done a small simulation study with R-code below. Feel free to play around with varying Ng and SSpg in the code, but the really important parameter is mTE, the mean effect size among those groups that do not correspond to true nulls.

With mTE <- 0, we find that indeed the true nulls yield significance at a nominal alpha-value of 0.05 less than 5% of the time. However, with mTE <- 5, and using the default values for Ng and SSpg in the code, these same groups produce Type 1 errors 100% of the time.

Code below deals with a random intercepts model for computational efficiency. The core conceptual issues are the same as in a random slopes model, I think.



type1 <- vector()
type2 <- vector()

Ng <- 400   # Must be an even number
            # We will split this number of groups in half.  One half will have true effect sizes of zero; these
              # are the true nulls.  The other half will have normally distributed effect sizes (these are
              # necessary to include in order to ensure that we estimate non-zero random effect variance).
SSpg = 20 # Sample size per group
mTE <- 0 # The mean of the true effect size for the groups that do not correspond to true nulls.

for(i in 1:100){
groupmeans <- c((rnorm(Ng/2)+mTE), rep(0,Ng/2))
# In this model, there are four hundred groups (to provide a reasonable sample size for calculating the FDR)
# The second 200 all have effect sizes of zero; these are the true nulls that might be subject to Type 1 error.
# The first 200 obey the random-effects specification. Including some groups like this is necessary to prevent
# the model from estimating zero random effect variance.

testdata <- as.data.frame(matrix(data = NA, nrow=Ng*SSpg, ncol=2))
colnames(testdata) <- c('group', 'y')

testdata$group <- rep(c(1:Ng), SSpg)
testdata$y <- rnorm(Ng*SSpg, groupmeans[testdata$group])

my.model <- lmer(y ~ (1|group), data = testdata)

# Get p-values for each group's effect using a normal approximation.

cV <- ranef(my.model, condVar = TRUE, drop=T)  
ranvar <- attr(cV[[1]], "postVar")
allvar <- ranvar + diag(vcov(my.model))
BLUPsd <- sqrt(allvar[1])

# check which group-level effects are significantly different from zero.
sig <- which(cV$group + fixef(my.model)  - 1.96*BLUPsd > 0 | cV$group + 1.96*BLUPsd < 0)
type1[i] <- length(which(sig > (Ng/2)))

# Expected number of false discoveries if p-values are uniformly distributed and independent over the 
# is Ng/40 true nulls

As Amoeba points out, the very idea of getting p-values from random effect BLUPs is not standard, as has been discussed here. The objections raised in the previous link are unfamiliar to me; from my (generally Bayesian) perspective I see nothing wrong with the idea of a credible interval around an individual group mean in a hierarchical model.

Bolker provides a sketch of how to get confidence intervals for the BLUPs from lme4, and the above computation of p-values is based on treating these "confidence intervals" uncritically as standard-fare confidence intervals. Briefly, we take the variance in the conditional mean for the BLUP (note that because this model is normal, the conditional mean and the conditional mode are one in the same), add to it the conditional variance in the intercept, and compute the z-score and p-value in just the way one might expect. I don't actually know whether the use of a normal distribution here is exact (because the whole model is Gaussian) or approximate. Note that some of the subtleties and pitfalls inherent in doing this actually involve covariances between random effects that are not present in this model. However, the approach used here does NOT account for the non-independence between the fixed intercept estimate and the random intercept estimates. I'm pretty sure that if anything, this would make my p-values anti-conservative, so the fact that the type 1 error rate is less than alpha (iff the true mean is close to zero) is still noteworthy.

  • 5
    $\begingroup$ I think the general philosophy of mixed models is that one does NOT perform inference on individual BLUPs of a random effect. If you are interested in individual values (in your case, group slopes), then use group as a fixed effect or simply conduct separate regressions for each group. Treating group as a random effect IMHO only makes sense when one is not interested in individual group values and considers between-group differences as a nuisance variability. $\endgroup$
    – amoeba
    Apr 27, 2018 at 21:49
  • $\begingroup$ In addition to @amoeba 's correct comment (+1), can you please set the seed of this code demonstration? In addition: I have the impression that 100 runs are not enough and probably a number of around 2000 would be more adequate. $\endgroup$
    – usεr11852
    Apr 27, 2018 at 22:48
  • $\begingroup$ @amoeba fair enough. I honestly don't know much about random effects models in a frequentist framework, but it's pretty common to hear Bayesians sing the praises of hierarchical models precisely to get better inference on the group-level effects by borrowing strength. Perhaps it would be more meaningful to ask the question purely from a Bayesian standpoint focused on how to interpret the credible intervals. $\endgroup$ Apr 28, 2018 at 1:21
  • $\begingroup$ @usεr11852 Note that each of the 100 runs contains 200 true null hypotheses embedded within it, so the simulation as a whole (as currently coded) computes p-values against 20,000 true null hypotheses. $\endgroup$ Apr 28, 2018 at 1:23
  • $\begingroup$ @JacobSocolar: Fair point, but we still get quite different end plot when running the same code twice. $\endgroup$
    – usεr11852
    Apr 28, 2018 at 11:42

1 Answer 1


If one of your model p-values has the meaning "probability of observing a parameter estimate that is as or more extreme than the one estimated, assuming the null hypothesis is true" (i.e. the standard interpretation of p-values), then you need to adjust for multiple comparisons.

This is because the meaning of "the probability of rejecting the null hypothesis when the null hypothesis is true" (given some preferred Type I error rate $\alpha$), is no longer coherent when there is no longer "the null hypothesis," but rather more than one null hypothesis.

That said: go with false discovery rate methods—almost certainly the Benjamini-Hochberg procedure (rather than the Benjamini-Yekutieli as a violation of positive dependence under regression is difficult to conceive)—as they (1) are not sensitive to/don't require a conceptually undefined "family", and (2) adapt to new information about the probability that a null hypothesis is true (i.e. if you have rejected some number of tests already, you should no longer believe 100% that all your remaining null hypotheses are true).

  • $\begingroup$ The question isn't "SHOULD one adjust for multiple comparisons", but rather "HOW should one adjust for multiple comparisons". Because we have a shrinkage estimator for the slopes, it seems possible at least in principle to control the FDR with a test that is uniformly more powerful than Benjamani-Hochberg. Based on the linked post from Gelman (linked in the OP), it appears that in some cases the necessary correction might be quite slight. $\endgroup$ Apr 26, 2018 at 4:02
  • $\begingroup$ @JacobSocolar Can you explain what a p-value's formal definition is in such models if it differs from the definition I gave? $\endgroup$
    – Alexis
    Apr 26, 2018 at 4:54
  • $\begingroup$ it does not differ from the definition you gave. $\endgroup$ Apr 26, 2018 at 6:05
  • $\begingroup$ @JacobSocolar Then my answer stands. However, I will edit it to include the Benjamini-Hochberg adjustment, explicitly. Let me know if you want to understand the steps required in B-H, and I will edit them in also. (And please have patience if I am still not understanding your question, I am receptive to my own learning :). $\endgroup$
    – Alexis
    Apr 26, 2018 at 15:01
  • 2
    $\begingroup$ I am afraid you have lost me. If you are not actually interested in inference using frequentist p-values (or frequentist CIs), then bets are off, and I am not sure what you are asking. multiple comparisons adjustments are specifically for thinking about frequentist inference given the meanings of $\alpha$, p and a Type I error. (I'm not qualified to comment on Bayesian approaches to Type I error in hypothesis testing) $\endgroup$
    – Alexis
    Apr 27, 2018 at 17:14

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.