# variance component possibly contributing to conflicting results in lme4 models

I have analyzed some data for a lab mate. The data has 4 levels (control, low, med, high), a random effect (rodent litter) and a continuous outcome (receptor binding in brain). He initially fit models with two levels, each of which was a comparison to the control in lme4:

mod_one <- lmer(OT_BNST ~ Tx + (1|Litter), data = males)
Random effects:
Groups   Name        Variance Std.Dev.
Litter   (Intercept) 33551    183.17
Residual              2159     46.47
Number of obs: 19, groups:  Litter, 14

Fixed effects:
Estimate Std. Error t value
(Intercept)   819.18      51.48  15.912
Tx2          -143.12      28.67  -4.992


This provided a significant effect. However, the estimates are very far from the empirical means:

control          medium
765.8853       736.9560


I argued that, based on only 14 "clusters" with only 19 data points, we are getting a very unstable estimate for the random intercept. In other words, 5 of the groups consist of n = 1 and these random intercepts are simply shrunken contrasts from the fixed effect intercept. Here are the conditional modes of the random intercepts:

$Litter (Intercept) 144-21 144.87999 144-23 -142.96796 150-21 -265.86454 151-11 103.74252 154-21 197.17572 154-22 -175.03226 163-11 63.97377 164-13 92.97553 165-2 -286.77684 166-2 136.25315 167-2 296.04793 168-2 10.16676 169-2 -131.51485 170-2 -43.05891  I then fit the model with all levels in the data (to me, more data can never be bad thing): Random effects: Groups Name Variance Std.Dev. Litter (Intercept) 2209 47 Residual 19325 139 Number of obs: 35, groups: Litter, 16 Fixed effects: Estimate Std. Error t value (Intercept) 768.35 44.14 17.407 Tx1 -55.23 63.70 -0.867 Tx2 -32.36 66.05 -0.490 Tx3 -12.08 68.60 -0.176  As can be seen the SD of Litter changes drastically and is now much more reasonable. Also, rather than 19 outcomes and 14 litters, there is now 35 outcomes from 16 litters. Moreover, the estimates are very close to the empirical means:  control low medium high 765.8853 706.9677 736.9560 755.8584  This model appears more reasonable especially when looking at the random intercepts: $Litter
(Intercept)
144-21  12.2291870
144-22 -17.3702645
144-23 -13.5126138
150-21 -30.6314011
151-11  23.7861574
154-21  23.2087960
154-22 -34.4604951
163-11   0.9603561
163-13   3.1420800
164-13  -8.0187069
165-2  -26.7697982
166-2   27.5962154
167-2   26.1703544
168-2    1.1041750
169-2    3.3305589
170-2    9.2353993


While I often use lme4, I have never experienced my estimates being affected so drastically as in mod_one. That is, the empirical difference is close to 30, while the model is estimating it at 150. However, I have also never really considered using a multilevel model with a very small sample as in mod_one and I generally use Bayesian methods (brms) where I can specify a prior on the SD of Litter.

Further complicating matters is that aov with an Error term and two treatment levels also provides a significant result:

aov(OT_BNST ~ Tx + Error(Litter), data = males)

Error: Litter
Df Sum Sq Mean Sq F value Pr(>F)
Tx         1  31315   31315   0.953  0.348
Residuals 12 394336   32861

Error: Within
Df Sum Sq Mean Sq F value  Pr(>F)
Tx         1  61575   61575   30.37 0.00529 **
Residuals  4   8110    2027


In contrast, here is the results from the model with all four levels:

 Error: Litter
Df Sum Sq Mean Sq F value Pr(>F)
Tx         3  44374   14791   0.519  0.677
Residuals 12 342085   28507

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Tx         3  16197    5399   0.306  0.821
Residuals 16 282578   17661


While I am confident with my non-significant model in lme4, I am not entirely sure what is going on other than, simply, estimating the random effect of Litter in mod_one is problematic. Moreover, calculating the sum of squares in aov for "blocking" variables, in this case Litter, is likely also problematic since 5 groups only consist of one and the sample sizes of Tx are unequal. Finally, I do not think it is reasonable to selectively do basically pairwise comparisons when everything can be specified in the same model.

I am hoping someone has some insight into what is going on and especially why the estimates of mod_one are so far from the empirical means.

• In addition, Cohen's d based on the raw data and for the two level model is - 0.17, [-1.12, 0.87]. – D_Williams Sep 4 '16 at 14:19
• Please post a reproducible example. Without that it is hard to say what might be going on. Certainly the variance of the random intercepts seems very high. Have you looked at their conditional modes with ranef() ? Perhaps this will give some insight. – Robert Long Sep 4 '16 at 15:32
• Robert: I agree about posting the data, but need to get the OK from lab mate as the data is not mine. Yes, the conditional modes are now in post. – D_Williams Sep 4 '16 at 15:40
• @RobertLong I have provided the random intercepts of Litter. – D_Williams Sep 4 '16 at 21:49