If the categorical variable is retained in my final model in R, then why does the post hoc analysis say the levels do not differ?

Question

I am performing model selection in R with the anova() function, and my categorical variable was maintained in my final model, but when I did a post hoc analysis with the emmeans() function, it told me the levels did not differ. What does it mean?

I use R software, and I am studying how the body condition of a species of fish varies in 3 kinds of rivers: preserved, slightly urban and very urban. Each category has one replicate, so that means I have 2 preserved rivers, 2 slightly urbanized rivers and 2 very urbanized rivers, which means that "river" is a random factor, and "category of urbanization" is my fixed factor and predictor variable with 3 levels. While performing model selection in R with the anova() function, the categorical variable "category" is maintained:

`#it is a linear mixed model because condition is normally distributed
> lmm.1 <- lmer(condition ~ category.of.urbanization + (1|river), data = fish) 
> lmm.null <- lmer(condition ~ 1 + (1|river), data = fish) 
> anova(lmm.null, lmm.1)
refitting model(s) with ML (instead of REML)
Data: fish
Models:
lmm.null: condition ~ 1 + (1 | river)
lmm.1: condition ~ category.of.urbanization + (1 | river)
            npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)   
lmm.null       3 -214.42 -205.37 110.21  -220.42                        
lmm.1          5 -219.80 -204.71 114.90  -229.80 9.3806  2   0.009184 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1`

` I had forgotten to include my model's summary in my my question, so here it is:

> summary(lmm.1)
Linear mixed model fit by REML ['lmerMod']
Formula: condition ~ category.of.urbanization + (1 | river)
Data: fish

REML criterion at convergence: -214.3

Scaled residuals: 
Min       1Q   Median       3Q      Max 
-2.65967 -0.54776 -0.07734  0.56748  2.79995 

Random effects:
Groups   Name        Variance Std.Dev.
river    (Intercept) 0.001965 0.04432 
Residual             0.012381 0.11127 
Number of obs: 151, groups:  river, 6

Fixed effects:
                                 Estimate Std. Error t value
(Intercept)                      -0.12808    0.04154  -3.084
category.of.urbanizationslightly  0.15972    0.05376   2.971
category.of.urbanizationvery      0.17063    0.05432   3.141

Correlation of Fixed Effects:
              (Intr) ctgr.d.rbnzcp
ctgr.d.rbnzcp -0.773              
ctgr.d.rbnzcm -0.765  0.591    

My p-value is 0.009184, meaning category of urbanization is an important predictor, and I expected that at least one level of the categorical variable would be different from the others. However, when trying to do a post hoc analysis, I called the emmeans() function, and R says that none of the levels differ, because the p-values are all above 0.05:

`> emmeans(lmm.1, pairwise ~ category.of.urbanization)
Registered S3 methods overwritten by 'broom':
method            from  
tidy.glht         jtools
tidy.summary.glht jtools
$emmeans
category.of.urbanization  emmean     SE   df lower.CL upper.CL
preserved                -0.1281 0.0441 3.42  -0.2592  0.00304
slightly urban            0.0316 0.0341 2.20  -0.1030  0.16632
very urban                0.0425 0.0350 2.43  -0.0852  0.17032

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
contrast                    estimate     SE   df t.ratio p.value
preserved - slightly urban   -0.1597 0.0558 2.83  -2.863  0.1324
preserved - very urban       -0.1706 0.0563 2.95  -3.028  0.1123
slightly urban - very urban  -0.0109 0.0489 2.32  -0.223  0.9733

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 3 estimates `

Please, what does this mean? How is the predictor variable significant, but with levels that aren't different? I have 151 fish, so my number of data and observations is not very low. I am sorry if I've made spelling mistakes, English is not my native language.

Answer 1

Some important information that you did not show is the between and within river variability (do a summary on the lmm.1 object). You have 151 fish, but only 6 rivers. If the between river variability is large relative to the within river variability then the 6 is the much more important number for the power of the inference.

Note also the degrees of freedom in the comparisons, they are very low (below 3). Remember that all the p-values are approximations once we go beyond simple models. You are computing approximate p-values based on the chi-squared distribution in one case and from a t distribution with fewer than 3 degrees of freedom in the other case. The square of a t is asymptotically chi-squared as the number of degrees of freedom go to infinity, so for large enough samples they should agree, but 3 d.f. is not large enough to expect agreement. I expect that the p-value based on the chi-squared is overly optimistic and those on the t a little bit pessimistic.

I would look at the estimates of the means as well as the differences and determine if those are potentially meaningful differences (practical significance), because if not, then there really is not reason to proceed. If they are potentially meaningful, then you either need more data (more rivers will be more important than more fish in the current rivers) or another way to bring in more information (possibly use a Bayesian analysis with informative priors).

Answer 2

Already without the mixed effects, you can get differences between different tests. Because

• The likelihood ratio test is approximated by a chi-squared distribution which is a bad approximation for small number of observations. You have 151 fish, but the degrees of freedom is effectively low because you only have three times two variations in the random effect.

• You have a situation where two values are very close to each other. In such case the p-value of an anova test will be relatively smaller (more significant) than the p-value of (corrected) individual t-tests.

  See:

  • Can ANOVA be significant when none of the pairwise t-tests is?
  • Discrepancy between the results of the ANOVA and the post-hoc test: How should be such results interpreted and presented?

Below is a simulation in R for a situation without mixed effects that shows how you can get the discrepancies.

You get different p-values

TEST                  p-value
Likelihood ratio test 0.011 
ANOVA test            0.103
TukeyHSD              0.129

Code:

### dummy data
d = 0.25
y = c(0,0,1,1,1,1) +
    c(-d,d,
      -d,d,
      -d,d)
x = as.factor(c(0,0,1,1,2,2))

mm = lm(y~x)

### LR test 0.01058
lmtest::lrtest(mm)

### ANOVA test 0.103
mod = aov(mm)
summary(mod)

### Tukey test 0.128822
TukeyHSD(mod)

Answer 3

There are two things going on:

1. You may have a lot of fish, but the experimental units for category.of.urbanization are the rivers, not the fish. For that reason, the degrees of freedom for each comparison is quite low. So although the $t$ ratios for comparing preserved with the other two conditions are similar to the model output, the degrees of freedom are both less than 3, and that makes the P values large.

2. In addition, as annotated in the emmeans output, the P values for the comparisons are Tukey-adjusted, and that inflates the P values, whereas the ones in the model summary are unadjusted (and also use the wrong d.f.)