2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ In addition, Cohen's d based on the raw data and for the two level model is - 0.17, [-1.12, 0.87]. $\endgroup$ – D_Williams Sep 4 '16 at 14:19
  • $\begingroup$ 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. $\endgroup$ – Robert Long Sep 4 '16 at 15:32
  • $\begingroup$ 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. $\endgroup$ – D_Williams Sep 4 '16 at 15:40
  • $\begingroup$ @RobertLong I have provided the random intercepts of Litter. $\endgroup$ – D_Williams Sep 4 '16 at 21:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.