Laurent Fanchon
2006-Apr-25 15:00 UTC
[R] lme: how to compare random effects in two subsets of data
Dear R-gurus, I have an interpretation problem regarding lme models. I am currently working on dog locomotion, particularly on some variation factors. I try to figure out which limb out of 2 generated more dispersed data. I record a value called Peak, around 20 times for each limb with a record. I repeat the records during a single day, and on several days. I tried to build two models, one for each limb : Dog.Left <- lme (fixed=Peak~1, data=Loco, subset=Limb=="Left",random=~1|Dog/Day/Record) Dog.Right <- lme (fixed=Peak~1, data=Loco, subset=Limb=="Right",random=~1|Dog/Day/Record) This allows to determine the variance attributable to each factor. Record represents the within-day variation, Day represents the between-day variation. This gives the following results : VarCorr (Dog.Left) Variance StdDev Dog = pdLogChol(1) (Intercept) 564.55587 23.760384 Day = pdLogChol(1) (Intercept) 54.63027 7.391229 Record = pdLogChol(1) (Intercept) 23.29377 4.826362 Residual 27.46464 5.240672 VarCorr(Dog.Right) Variance StdDev Dog = pdLogChol(1) (Intercept) 552.11246 23.497074 Day = pdLogChol(1) (Intercept) 70.72088 8.409571 Record = pdLogChol(1) (Intercept) 21.94594 4.684649 Residual 29.68476 5.448373 This shows that the variance might be different for each limb. For example, the variance attributable to Day might be higher for the Right limb. This is the first part of my interpretation, and I hope to be right. What do you think?? Then, the question is : are these differences statistically significant. I am not sure of how to investigate this question. I tried to compare several models : model1<- lme (fixed=Peak~Limb, data=Loco, random=list(Dog=~Limb,Day=~Limb,Record=~Limb)) this is the more complicated model model2<-lme (fixed=Peak~Limb, data=Loco, random=list(Dog=~Limb,Day=~Limb,Record=~1)) anova (model1,model2) showed no difference model3<-lme (fixed=Peak~Limb, data=Loco, random=list(Dog=~Limb,Day=~1,Record=~1)) anova (model2,model3) showed a significant difference <0.0001 model2 seems to be the best model. Does it means the difference of variance between the two limb is significant for between-day variation and is unsignificant for within-day variation?? Finally VarCorr (model2) gives : Variance StdDev Corr Dog = pdLogChol(Limb) (Intercept) 567.553021 23.823371 (Intr) LimbRight 7.249064 2.692409 -0.166 Day = pdLogChol(Limb) (Intercept) 53.888346 7.340868 (Intr) LimbRight 4.863394 2.205310 0.363 Record = pdLogChol(1) (Intercept) 22.418031 4.734768 Residual 28.740451 5.361012 I am not sure to understand this issue. The global variance attributable to Day is roughly 53.88 (random effect on the intercept). And the differences between the two limbs might be increased according to a variance of 4.86 (random effect on the slope). Is that right? But does this also make it possible to determine which limb had the highest variance? I guess if I change the order of the Limb factor (Right<Left) I will get the same results with LimbLeft. How to determine the limb with the highest variance? Do I have to refer to the first two models (Variance attributable to Day was a little higher for the Right) ? Any help will be very appreciated, Laurent Fanchon Ecole Nationale V?t?rinaire d'Alfort National Vet School of Alfort (France)
Spencer Graves
2006-May-02 16:14 UTC
[R] lme: how to compare random effects in two subsets of data
(comments in line) Laurent Fanchon wrote:> Dear R-gurus, > > I have an interpretation problem regarding lme models. > > I am currently working on dog locomotion, particularly on some variation > factors. > I try to figure out which limb out of 2 generated more dispersed data. > > I record a value called Peak, around 20 times for each limb with a record. > I repeat the records during a single day, and on several days. > > I tried to build two models, one for each limb : > Dog.Left <- lme (fixed=Peak~1, data=Loco, > subset=Limb=="Left",random=~1|Dog/Day/Record) > Dog.Right <- lme (fixed=Peak~1, data=Loco, > subset=Limb=="Right",random=~1|Dog/Day/Record) > > This allows to determine the variance attributable to each factor. > Record represents the within-day variation, Day represents the > between-day variation. > > This gives the following results : > VarCorr (Dog.Left) > Variance StdDev > Dog = pdLogChol(1) > (Intercept) 564.55587 23.760384 > Day = pdLogChol(1) > (Intercept) 54.63027 7.391229 > Record = pdLogChol(1) > (Intercept) 23.29377 4.826362 > Residual 27.46464 5.240672 > > VarCorr(Dog.Right) > Variance StdDev > Dog = pdLogChol(1) > (Intercept) 552.11246 23.497074 > Day = pdLogChol(1) > (Intercept) 70.72088 8.409571 > Record = pdLogChol(1) > (Intercept) 21.94594 4.684649 > Residual 29.68476 5.448373 > > This shows that the variance might be different for each limb. > For example, the variance attributable to Day might be higher for the > Right limb. > > This is the first part of my interpretation, and I hope to be right. > What do you think?? > > Then, the question is : are these differences statistically significant. > I am not sure of how to investigate this question. > > I tried to compare several models : > model1<- lme (fixed=Peak~Limb, data=Loco, > random=list(Dog=~Limb,Day=~Limb,Record=~Limb)) this is the more > complicated model > model2<-lme (fixed=Peak~Limb, data=Loco, > random=list(Dog=~Limb,Day=~Limb,Record=~1)) > anova (model1,model2) showed no difference > model3<-lme (fixed=Peak~Limb, data=Loco, > random=list(Dog=~Limb,Day=~1,Record=~1)) > anova (model2,model3) showed a significant difference <0.0001 > > model2 seems to be the best model. > Does it means the difference of variance between the two limb is > significant for between-day variation and is unsignificant for > within-day variation??NO, it does NOT mean that the variance between the two limbs is significantly different between days. The model includes additive errors for each day as e[day]+e[day & limb], and the contribution from e[day & limb] is statistically significant.> > Finally VarCorr (model2) gives : > Variance StdDev Corr > Dog = pdLogChol(Limb) > (Intercept) 567.553021 23.823371 (Intr) > LimbRight 7.249064 2.692409 -0.166 > Day = pdLogChol(Limb) > (Intercept) 53.888346 7.340868 (Intr) > LimbRight 4.863394 2.205310 0.363This estimates var(e[day]) = 53.9 and var(e[day & limb]) = 4.9. To test a difference in variance between limbs, please see Sec. 5.2 in Pinheiro and Bates (2000) Mixed-Effects Models in S and S-PLUS (Springer). I couldn't find the electronic catalog on the web site for the Biblioth?que de l'?cole nationale v?t?rinaire d?Alfort, so I couldn't check to see if they have a copy. If they don't, please feel free to advise them for me that you've heard that this is the best book on this subject available today, and it should find many users among researchers such as yourself. In my opinion, Doug Bates is the leading expert on variance components questions in the world today, and an institution such as the ?cole nationale v?t?rinaire d?Alfort should have a copy of this book in their library, and their researchers who routinely must work with these kinds of data should make routine use of this book, if they don't already. hope this helps, spencer graves> Record = pdLogChol(1) > (Intercept) 22.418031 4.734768 > Residual 28.740451 5.361012 > > I am not sure to understand this issue. > The global variance attributable to Day is roughly 53.88 (random effect > on the intercept). And the differences between the two limbs might be > increased according to a variance of 4.86 (random effect on the slope). > Is that right? > But does this also make it possible to determine which limb had the > highest variance? I guess if I change the order of the Limb factor > (Right<Left) I will get the same results with LimbLeft. How to determine > the limb with the highest variance? > Do I have to refer to the first two models (Variance attributable to Day > was a little higher for the Right) ? > > Any help will be very appreciated, > > Laurent Fanchon > Ecole Nationale V?t?rinaire d'Alfort > National Vet School of Alfort (France) > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html