[R] FW: Re: Doubt about nested aov output
Your response nicely clarifies a question that I've had for a long time, but which I've dealt with by giving each subject a unique label. Unless I'm missing something, both techniques should work as the toy example below gives exactly the same output in all 3 cases below (forgetting about the convergence problem). Would there be a reason to prefer labeling the levels one way or another or is it just a matter of convenience? library(lmer) y - rnorm(15) cond - gl(3, 5, 15) obs - gl(15, 1) subj - gl(5, 1, 15) dd - data.frame(y = y, cond = cond, obs = obs, subj = subj) l1 - lmer(y~cond + (1|cond:obs), dd) l2 - lmer(y~cond + (1|cond:subj), dd) l3 - lmer(y~cond + (1|obs), dd) Douglas Bates a écrit: The difference between models like lmer(Glycogen~Treatment+(1|Rat)+(1|Rat:Liver)) and lmer(Glycogen~Treatment+(1|Treatment:Rat)+(1|Treatment:Rat:Liver)) is more about the meaning of the levels of Rat than about the meaning of Treatment. As I understood it there are three different rats labelled 1. There is a rat 1 on treatment 1 and a rat 1 on treatment 2 and a rat 1 on treatment 3. Thus the levels of Rat do not designate the experimental unit, it is the levels of Treatment:Rat that do this. -- Ken Knoblauch Inserm U371 Cerveau et Vision Dept. of Cognitive Neuroscience 18 avenue du Doyen Lépine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: +33 (0)6 84 10 64 10 http://www.lyon.inserm.fr/371/ __ R-help@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
Re: [R] FW: Re: Doubt about nested aov output
On 9/8/05, Ken Knoblauch [EMAIL PROTECTED] wrote: Your response nicely clarifies a question that I've had for a long time, but which I've dealt with by giving each subject a unique label. Unless I'm missing something, both techniques should work as the toy example below gives exactly the same output in all 3 cases below (forgetting about the convergence problem). Would there be a reason to prefer labeling the levels one way or another or is it just a matter of convenience? library(lmer) y - rnorm(15) cond - gl(3, 5, 15) obs - gl(15, 1) subj - gl(5, 1, 15) dd - data.frame(y = y, cond = cond, obs = obs, subj = subj) l1 - lmer(y~cond + (1|cond:obs), dd) l2 - lmer(y~cond + (1|cond:subj), dd) l3 - lmer(y~cond + (1|obs), dd) I prefer to have a grouping factor constructed with unique levels for each distinct unit. The only reason I mention constructions like Treatment:Rat in the original part of this thread is that data are often provided in that form. Reusing subject labels within another group is awkward and can be error prone. One of the data sets I examine in the MlmSoftRev vignette of the mlmRev package is called Exam and has student identifiers within schools. The student identifiers are not unique but the school:student combination should be. It isn't. These data have been analyzed in scores of books and articles and apparently none of the other authors bothered to check this. There are some interesting ramifications such as some of the schools that are classified as mixed-sex schools are likely single-sex schools because the only student of one of the sexes in that school is apparently mislabelled. BTW, in your example you have only one observation per level of 'obs' so you can't use obs as a grouping factor as this variance component would be completely confounded with the per-observation noise. Douglas Bates a écrit: The difference between models like lmer(Glycogen~Treatment+(1|Rat)+(1|Rat:Liver)) and lmer(Glycogen~Treatment+(1|Treatment:Rat)+(1|Treatment:Rat:Liver)) is more about the meaning of the levels of Rat than about the meaning of Treatment. As I understood it there are three different rats labelled 1. There is a rat 1 on treatment 1 and a rat 1 on treatment 2 and a rat 1 on treatment 3. Thus the levels of Rat do not designate the experimental unit, it is the levels of Treatment:Rat that do this. -- Ken Knoblauch Inserm U371 Cerveau et Vision Dept. of Cognitive Neuroscience 18 avenue du Doyen Lépine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: +33 (0)6 84 10 64 10 http://www.lyon.inserm.fr/371/ __ R-help@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
Re: [R] FW: Re: Doubt about nested aov output
Thank you for your response. The single response/observer most probably explains the complaints that lmer was giving for my example. Maybe this small modification provides a better example and corrects a more serious error in my previous post: library(lme4) y-rnorm(30) cond - rep(gl(3,5,15), 2) obs-rep(gl(15,1), 2) subj-rep(gl(5,1,15), 2) dd-data.frame(y=y,cond=cond,obs=obs,subj=subj) l1 - lmer(y~cond + (1|cond:obs), data=dd) l2 - lmer(y~cond + (1|cond:subj), data=dd) l3 - lmer(y~cond + (1|obs), dd) Understanding the notation is often about 99% of the job, and it is very helpful to see multiple ways to accomplish the same thing. Douglas Bates a écrit: I prefer to have a grouping factor constructed with unique levels for each distinct unit. The only reason I mention constructions like Treatment:Rat in the original part of this thread is that data are often provided in that form. Reusing subject labels within another group is awkward and can be error prone. One of the data sets I examine in the MlmSoftRev vignette of the mlmRev package is called Exam and has student identifiers within schools. The student identifiers are not unique but the school:student combination should be. It isn't. These data have been analyzed in scores of books and articles and apparently none of the other authors bothered to check this. There are some interesting ramifications such as some of the schools that are classified as mixed-sex schools are likely single-sex schools because the only student of one of the sexes in that school is apparently mislabelled. BTW, in your example you have only one observation per level of 'obs' so you can't use obs as a grouping factor as this variance component would be completely confounded with the per-observation noise. Douglas Bates a écrit: The difference between models like lmer(Glycogen~Treatment+(1|Rat)+(1|Rat:Liver)) and lmer(Glycogen~Treatment+(1|Treatment:Rat)+(1|Treatment:Rat:Liver)) is more about the meaning of the levels of Rat than about the meaning of Treatment. As I understood it there are three different rats labelled 1. There is a rat 1 on treatment 1 and a rat 1 on treatment 2 and a rat 1 on treatment 3. Thus the levels of Rat do not designate the experimental unit, it is the levels of Treatment:Rat that do this. -- Ken Knoblauch Inserm U371 Cerveau et Vision Dept. of Cognitive Neuroscience 18 avenue du Doyen Lépine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: +33 (0)6 84 10 64 10 http://www.lyon.inserm.fr/371/ -- Ken Knoblauch Inserm U371 Cerveau et Vision Dept. of Cognitive Neuroscience 18 avenue du Doyen Lépine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: +33 (0)6 84 10 64 10 http://www.lyon.inserm.fr/371/ __ R-help@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
[R] FW: Re: Doubt about nested aov output
Ronaldo, Further to my previous posting on your Glycogen nested aov model. Having read Douglas Bates' response and Reflected on his lmer analysis output of your aov nested model example as given.The Glycogen treatment has to be a Fixed Effect.If a 'treatment' isn't a Fixed Effect what is ? If Douglas Bates' lmer model is modified to treat Glycogen Treatment as a purely Fixed Effect,with Rat and the interaction Rat:Liver as random effects then-- model.lmer-lmer(Glycogen~Treatment+(1|Rat)+(1|Rat:Liver)) summary(model.lmer) Linear mixed-effects model fit by REML Formula: Glycogen ~ Treatment + (1 | Rat) + (1 | Rat:Liver) AIC BIClogLik MLdeviance REMLdeviance 239.095 248.5961 -113.5475 238.5439 227.095 Random effects: GroupsNameVariance Std.Dev. Rat:Liver (Intercept) 2.1238e-08 0.00014573 Rat (Intercept) 2.0609e+01 4.53976242 Residual 4.2476e+01 6.51733769 # of obs: 36, groups: Rat:Liver, 6; Rat, 2 Fixed effects: Estimate Std. Error DF t value Pr(|t|) (Intercept) 140.5000 3.7208 33 37.7607 2.2e-16 *** Treatment2 10.5000 2.6607 33 3.9463 0.0003917 *** Treatment3 -5. 2.6607 33 -2.0045 0.0532798 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) Trtmn2 Treatment2 -0.358 Treatment3 -0.358 0.500 anova(model.lmer) Analysis of Variance Table Df Sum Sq Mean Sq Denom F valuePr(F) Treatment 2 1557.56 778.78 33.00 18.335 4.419e-06 *** which agrees with the aov model below. model - aov(Glycogen~Treatment+Error(Rat/Liver)) summary(model) John -Original Message- From: John Wilkinson (pipex) [mailto:[EMAIL PROTECTED] Sent: 07 September 2005 12:04 PM To: Ronaldo Reis-Jr. Cc: r-help Subject: Re: [R] Doubt about nested aov output Ronaldo , It looks as though you have specified your model incorrectly. In the Rats example ,the Treatment is the only fixed effect,Rat and Liver are random effects In aov testing for sig of 'Means' of Random Effects is pointless and that is why 'p' values are not given.Further more the interaction between a Random Effect and a Fixed Effect is also a Random Effect. The 'aov' with error structure terms output reflects this by only giving 'p' values to Fixed Effects and their interactions model - aov(Glycogen~Treatment+Error(Rat/Liver)) summary(model) Error: Rat Df Sum Sq Mean Sq F value Pr(F) #Rat is random effect Residuals 1 413.44 413.44 Error: Rat:Liver #Rat:Liver is Random effect Df Sum Sq Mean Sq F value Pr(F) Residuals 4 164.444 41.111 Error: Within Df Sum Sq Mean Sq F valuePr(F) Treatment 2 1557.56 778.78 18.251 8.437e-06 *** #Fixed effect Residuals 28 1194.78 42.67 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 I hope that this is of help. John __ R-help@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
