[R-sig-ME] Random effects correlations and factor reference levels

Wed Oct 14 12:17:16 CEST 2015

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

> On 13 Oct 2015, at 16:08, paul debes <paul.debes at utu.fi> wrote:
> 
> 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:
> 
>> 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)
>> 
>>> M1 <- lme(fixed = score ~ Machine, random = ~Machine|Worker, data=Machines)
>>> VarCorr(M1)
>> 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:
>> 
>>> Machines$Machine <- relevel(Machines$Machine, ref = "B")
>>> M2 <- update(M1)
>>> VarCorr(M2)
>> 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
>> ------------------------------------------------------------------------
>> 
>> 
>> 
>> 
>> 
>> 	[[alternative HTML version deleted]]
>> 
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 
> -- 
> Paul Debes
> DFG Research Fellow
> University of Turku
> Department of Biology
> Itäinen Pitkäkatu 4
> 20520 Turku
> Finland
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models