Skip to content

Random effects correlations and factor reference levels

3 messages · paul debes, Christopher Knight

Original
Christopher Knight
# 
I am fitting mixed effects models where there is a random effect which acts differently on each level of a factor, much as in the Machines example from the nlme package:

library(nlme)
Worker = pdLogChol(Machine)
            Variance   StdDev    Corr
(Intercept) 16.6405306 4.0792806 (Intr) MachnB
MachineB    34.5466908 5.8776433  0.484
MachineC    13.6150244 3.6898543 -0.365  0.297
Residual     0.9246296 0.9615766

I was interpreting these correlations to mean that the estimate of workers? scores on Machine B (or more precisely their BLUPS, interpreted as something like their residuals from the average on Machine B) are weakly positively correlated with their scores on Machine C (r=0.297). Because Machine A is the reference level of the factor ?Machine? (with treatment contrasts), I was also interpreting the ?(Intercept)? to represent Machine A, i.e. that the scores of the workers on Machine A are somewhat positively correlated with their scores on Machine B (r=0.484), but negatively with their scores on Machine C (r=-0.365).

If this were true, I would expect these conclusions to be robust to exactly which machine is the reference level in the contrasts for the Machine factor, but apparently they are not:
Worker = pdLogChol(Machine)
            Variance   StdDev    Corr
(Intercept) 74.3956398 8.6252907 (Intr) MachnA
MachineA    34.5466902 5.8776433 -0.910
MachineC    35.2950235 5.9409615 -0.882  0.805
Residual     0.9246296 0.9615766

What is going on here/what am I misinterpreting?

As expected, these two models appear to be equivalent by anova(M1,M2), yet the correlation among workers between machines A and B appears to be strongly negative in the second case (r = -0.91), where it seems moderately positive in the first case (r=0.484) and intervals() on the two models suggests no overlap of (admittedly pretty broad) confidence intervals in this case. Add to that the differences in the estimated variance components/standard deviations ? each model has an identical residual, but very different individual values (that for Machine C going from 13.6 to 35) and totals; at the same time the estimate for the variance among workers on Machine B in model M1 (34.5466908)  is suspiciously similar (albeit not quite identical) to the estimate for variance among workers on Machine A in model M2 (34.5466902).

The values are pretty much the same fitting with lmer (though in practice I need the correlation and weights arguments available with lme), which makes it feel like a problem with my interpretation. Or is it really just a case of too little data to resolve such ?subtle parameters clearly?

Dropping the fixed effect doesn?t make things much better and, on the models I am running in anger, there is more data and the differences are even more pronounced. It looks a bit like this issue: http://stats.stackexchange.com/questions/82102/changing-reference-level-for-contrasts-changes-results-in-r-3-0-2-lme4-1-1-2-vs
which was framed in terms of versions but, given that I get very similar issues using nlme and lme4, it doesn?t seem likely to be a version issue here (albeit I don?t have old versions immediately to hand to test this on).

Any insights much appreciated,

Chris

running nlme_3.1-122 in R 3.2.2 on OS X 10.10.5


------------------------------------------------------------------------
Dr Christopher Knight                             Michael Smith Building
Lecturer                                        Faculty of Life Sciences
Tel:  +44 (0)161 2755378                    The University of Manchester
room B.2012                                                  Oxford Road
tinyurl.com/knightFLS/<http://tinyurl.com/knightFLS/>                                Manchester M13 9PT
? . ,,><(((?>                                                         UK
------------------------------------------------------------------------
paul debes
# 
Hi Chris,

do you really need the treatment contrast?
I think the results based on the treatment contrasts for the random  
effects are more difficult to interpret than what you suggested (I can't  
interpret them either). If you simply remove the intercept: "random =  
~Machine-1|Worker" your interpretations of the worker-score correlations  
between machine levels hold and you get the actual variance estimates for  
each level of machine.

Best,
Paul




On Tue, 13 Oct 2015 16:57:32 +0300, Christopher Knight
<Chris.Knight at manchester.ac.uk> wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Christopher Knight
    

    

      
      #
    

  


  

      
Thank you, that does indeed make a whole lot more sense!

It would be good to get my head around what the intercept was doing originally (certainly not adding any extra parameters, which is perhaps why I hadn?t tried it), but it seems that that solves the interpretation problem in my real data too (albeit that, unlike the simpler example, the precise numbers aren?t entirely stable to re-specifying the reference level, but changes are associated with changes in the log likelihood and can be tweaked by changing the optimiser).

Thanks again,

Chris