Using 95% confidence intervals for pairwise comparisons in mixed effects model

Question

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) 
newdat 

Here were the results:

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

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

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

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.