R help - Dec 2007 - hierarchical linear models, mixed models and lme

Dear R-users,

I am trying to analyse the data of the box 10.5 in the Biometry from 
Sokal and Rohlf (2001) using R. This is a three-level nested anova with 
equal sample size : 3 different treatments are compared ; 2 rats (coded 
1 or 2) / treatment are studied ; 3 preparations (coded 1, 2 or 3) / 
rats are available ; 2 readings of the glycogen content  / preparations 
are realised. Treatment is fixed whereas Rats (nested in Treatment)  and 
Prep (nested in Rats) are random effects.

According to a previous discussion found in the R-help archives (January 
2007), I have tried the following formula :
>  box105.lme<-lme(content~treatment, box105.gd, random=~1|rats/prep)
However,  the formula summary(box105.lme) gives wrong estimates for the 
variance components (rats within treatments, prep within rats, readings 
within preps) ! Moreover the numbers of rats and preps are also wrong, 
with respectively 2 and 6 instead of 3*2=6 and 6* 3=18 !

I have also tried the following formula :
>  box105bis.lme<-lme(content~treatment, box105.gd, 
random=~1|treatment/rats/prep
In this case, the variance components as well as the number of rats and 
preps are correct. Nevertheless, I have two new problems : (1) the 
treatment is first treated as a random effects although it is fixed ! 
(2) there is a serious problem of df when treatment is then treated as a 
fixed effect (18 df for the intercept and 0 for the two other treatments !)

What's wrong ? I didn't find such design and data in Pinheiro and Bates 
(2000)

Thanks for your help,

Nicolas

Nicolas Ris <Nicolas.Ris <at> sophia.inra.fr> writes:
> I am trying to analyse the data of the box 10.5 in the Biometry from 
> Sokal and Rohlf (2001) using R. This is a three-level nested anova with 
> equal sample size : 3 different treatments are compared ; 2 rats (coded 
> 1 or 2) / treatment are studied ; 3 preparations (coded 1, 2 or 3) / 
> rats are available ; 2 readings of the glycogen content  / preparations 
> are realised. Treatment is fixed whereas Rats (nested in Treatment)  and 
> Prep (nested in Rats) are random effects.
> 
> According to a previous discussion found in the R-help archives (January 
> 2007), I have tried the following formula :
> >  box105.lme<-lme(content~treatment, box105.gd, random=~1|rats/prep)
> However,  the formula summary(box105.lme) gives wrong estimates for the 
Since your factors are codes as numbers, the first thing I would check if these
are really factors in your data frame. If these are numeric, you will get a
regression and degrees of freedom are totally off. Try

box105.gd$rats = as.factor(box105.gd$rats)
box105.gd$prep = as.factor(box105.gd$prep)

I am sure there will be some differences to Sokal/Rohlf afterwards, but you
should come closer. Please post a complete, self-running examples if you have
more questions.

Dieter