I'm running mixed effects models in R using Treatment and Species as crossed predictors and Group as a random effect using the nlme package, like so:

model <- lme(growth_rate ~ Treatment * Species, random = ~1|Group, data = data, method="REML")

I wanted to run a post-hoc test to perform pairwise comparisons, i.e. compare how different treatment/species combinations compare to each other.

I was having some trouble getting the emmeans package to run (though I have solved this issue now). Someone suggested an alternative to emmeans, which was to construct 95% confidence intervals for my predicted values (as described here). I was told that treatment/species combinations with confidence intervals that do not overlap could be considered to be significantly different. Is this correct and if so, why don't we need to correct for multiple comparisons?

Here is how I manually computed the confidence intervals from the model based on the code from the link:

newdat <- expand.grid(Treatment=c("A", "B", "C", "C", "E"), Species=c("X", "Y")) 
newdat$pred <- predict(model, newdat, level = 0)
## [-2] drops response from formula
Designmat <- model.matrix(formula(model)[-2], newdat)
##compute XVX' to get variance-covariance matrix of predictions and extract the diagonal to get variances of predictions
predvar <- diag(Designmat %*% vcov(model) %*% t(Designmat)) 
##take square root of variances to get the standard deviations (errors) of the predictions
newdat$SE <- sqrt(predvar)
cmult <- 1.96
newdat<- newdat %>% 
      mutate(CI_low = pred-cmult*SE) %>% 
      mutate(CI_high = pred+cmult*SE) 

Here were the results: Results from manually computed CIs

Then I tried running emmeans on the model to compare my results:

    t <- emmeans(model, pairwise ~ Treatment * Species, adjust = "tukey")

emmeans CIs

The results were different enough (no significant pairwise comparisons) that it changes my conclusions substantially.

  • 1
    $\begingroup$ In this case you do not need to correct for multiple comparisons. The levels of the factors (species and treatment) are accounted for by lme4 in the mixed-effects model model fitting, and you are using the model to predict. Different prediction methods may return slightly different estimates of predicted means and confidence intervals, but generally if intervals do not overlap then differences are statistically significant. Check the ggpredict and sjplot packages for some neat tools that can help compute these and visualize them. ggpredict::ggmeans seems appropriate here, which uses emmeans. $\endgroup$
    – hoganhaben
    Commented Nov 28, 2023 at 18:47
  • $\begingroup$ Ok, thank you so much! I guess I still don't quite understand what is going on under the hood- how the model is correcting for the levels of the factors? I ran emmeans with Tukey's because I was curious, and it appears to be much more conservative than using the confidence intervals the way I computed them manually. Why would that be the case? Thanks for the advice on ggpredict and sjplot- I'll look into those! Also, sorry I had the link slightly wrong- it is here: bbolker.github.io/mixedmodels-misc/… $\endgroup$
    – mels
    Commented Nov 28, 2023 at 19:42
  • 3
    $\begingroup$ Non-overlap of 95% CI is not the same thing as a significant difference at p < 0.05; see here for a discussion in the context of a t-test. The emmeans package is well vetted, so if your CI comparison is substantially different from what emmeans is telling you I fear that there might be a problem. Consider editing this question to show both types of analyses. The emmeans calculations only use the fixed-effect coefficients, but that's usually what you want to evaluate once you've handled Group with random intercepts. $\endgroup$
    – EdM
    Commented Nov 28, 2023 at 19:53
  • 2
    $\begingroup$ @hoganhaben Unfortunately, this case is not "in general," because confidence intervals are not independent in this circumstance. The simple approximations and considerations involved in comparing independent confidence intervals don't apply. $\endgroup$
    – whuber
    Commented Nov 28, 2023 at 19:55
  • 2
    $\begingroup$ The emmeans confidence intervals are t-based, not Z. Use qt(.975, 9) $\approx$ 2.262 (for df=9) or 2.571 (for df=5) where you used 1.96 (cmult) and they should match. I do want to reiterate EdM's point that the confidence intervals of the marginal means do not say much about the difference between them - non-overlap is a much more stringent requirement than the CI of their difference excluding zero. $\endgroup$
    – PBulls
    Commented Nov 28, 2023 at 21:06

1 Answer 1


The point estimates and standard errors agree perfectly between the two methods, and the confidence limits calculated the two ways are very similar (probably due to the difference between z- based and t-based intervals, as in a comment from @PBulls).

The correction for multiple comparisons is the real difference here. You need to make such corrections, as explained in an emmeans vignette. The temptation to find some point estimates that are outside of the confidence interval for some other estimate can be strong, but it must be avoided. Just looking for a few "big differences" among many possible differences is data dredging that will typically lead you astray.

Proper statistical tests on differences between estimates need to start with the estimates and the coefficient covariance matrix, applying the formula for the variance of a weighted sum of correlated variables. The emmeans package can do all that for you. It's certainly instructive to try constructing some tests of comparisons by hand, but if you do them correctly you shouldn't get different results from corresponding emmeans functions--for the standard errors of the estimated differences.

Then, however, you need to correct for multiple comparisons. With 10 estimated marginal means, you have 45 pairwise comparisons. Each additional comparison makes it harder to detect a true difference after proper correction.

Do you really care about all 45 comparisons?

Maybe all you care about are comparisons of each of 4 treatments against some standard treatment (maybe Treatment A?), within each species. That's only 8 comparisons.

Maybe you are willing to settle for a report that combines results for the two species. That ignores the interaction term and can be misleading in general; if you try that with emmeans you will probably get a corresponding warning along with the calculations. In your case, however, for each treatment you have the same number of observations for both species, so that might be OK. You would then have only 5 estimated marginal means, for only 10 all-pairwise comparisons and only 4 comparisons of each of 4 treatments against a standard treatment.

  • $\begingroup$ Thank you, this is really helpful. I can settle for comparisons within species. $\endgroup$
    – mels
    Commented Nov 29, 2023 at 20:39
  • $\begingroup$ Just to confirm, is this the correct code for comparing treatments only within a species: emmeans(model, specs = pairwise ~ Treatment|Species, adjust = "tukey") ? $\endgroup$
    – mels
    Commented Nov 30, 2023 at 6:57
  • 1
    $\begingroup$ @mels as I recall, that should give all 10 pairwise comparisons among the 5 treatments with Tukey corrections, separately within each species. I might be wrong; I often takes me a few trials to get things right. If you only want comparisons of multiple treatments individually against a standard treatment, use a form of trt.vs.ctrl instead. Note the caution against the pairwise~ syntax in the FAQ; it's safest to get the grid for the modeled data first and then do contrasts based on the grid. $\endgroup$
    – EdM
    Commented Nov 30, 2023 at 15:49

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.