[R] Doubt about nested aov output

[Douglas Bates]

On 9/6/05, Ronaldo Reis-Jr. <chrysopa at gmail.com> wrote:
> Hi Spencer,
> 
> Em Dom 04 Set 2005 20:31, Spencer Graves escreveu:
> >         Others may know the answer to your question, but I don't.  However,
> > since I have not seen a reply, I will offer a few comments:
> >
> >         1.  What version of R are you using?  I just tried superficially
> > similar things with the examples in ?aov in R 2.1.1 patched and
> > consistently got F and p values.
> 
> I'm using the R version 2.1.1 on Linux Debian
> Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0
> 
> >         2.  My preference for this kind of thing is to use lme in
> > library(nlme) or lmer in library(lme4).  Also, I highly recommend
> > Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer).
> 
> Yes, this is my preference too, but I need aov for classes.
> 
> >         3.  If still want to use aov and are getting this problem in R 2.1.1,
> > could you please provide this list with a small, self contained example
> > that displays the symptoms that concern you?  And PLEASE do read the
> > posting guide! "http://www.R-project.org/posting-guide.html".  It might
> > increase the speed and utility of replies.
> >
> >         spencer graves
> 
> I send the complete example. This is a example from the Crwaley's book
> (Statistical Computing: An introdution to data analysis using S-Plus.
> 
> This is a classical experiment to show pseudoreplication, from Sokal and Rohlf
> (1995).
> 
> In this experiments, It have 3 treatmens applied to 6 rats, for each rat it
> make 3 liver preparation and for each liver it make 2 readings of glycogen.
> This generated 6 pseudoreplication per rat. I'm interested on the effect os
> treatment on the glycogen readings.
> 
> Look the R analyses:
> 
> --------------------
> > Glycogen <-
> c(131,130,131,125,136,142,150,148,140,143,160,150,157,145,154,142,147,153,151,155,147,147,162,152,134,125,138,138,135,136,138,140,139,138,134,127)
> > Glycogen
>  [1] 131 130 131 125 136 142 150 148 140 143 160 150 157 145 154 142 147 153
> 151
> [20] 155 147 147 162 152 134 125 138 138 135 136 138 140 139 138 134 127
> > Treatment <- factor(rep(c(1,2,3),c(12,12,12)))
> > Treatment
>  [1] 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3
> Levels: 1 2 3
> > Rat <- factor(rep(rep(c(1,2),c(6,6)),3))
> > Rat
>  [1] 1 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 2
> Levels: 1 2
> > Liver <- factor(rep(rep(c(1,2,3),c(2,2,2)),6))
> > Liver
>  [1] 1 1 2 2 3 3 1 1 2 2 3 3 1 1 2 2 3 3 1 1 2 2 3 3 1 1 2 2 3 3 1 1 2 2 3 3
> Levels: 1 2 3
> >
> > ### Model made identical to the book
> >
> > model <- aov(Glycogen~Treatment/Rat/Liver+Error(Treatment/Rat/Liver))
> >
> > summary(model)
> 
> Error: Treatment
>           Df  Sum Sq Mean Sq
> Treatment  2 1557.56  778.78
> 
> Error: Treatment:Rat
>               Df Sum Sq Mean Sq
> Treatment:Rat  3 797.67  265.89
> 
> Error: Treatment:Rat:Liver
>                     Df Sum Sq Mean Sq
> Treatment:Rat:Liver 12  594.0    49.5
> 
> Error: Within
>           Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 18 381.00   21.17
> >
> > ### Model made by myself, I'm interested only in Treatment effects
> >
> > model <- aov(Glycogen~Treatment+Error(Treatment/Rat/Liver))
> >
> > summary(model)
> 
> Error: Treatment
>           Df  Sum Sq Mean Sq
> Treatment  2 1557.56  778.78
> 
> Error: Treatment:Rat
>           Df Sum Sq Mean Sq F value Pr(>F)
> Residuals  3 797.67  265.89
> 
> Error: Treatment:Rat:Liver
>           Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 12  594.0    49.5
> 
> Error: Within
>           Df Sum Sq Mean Sq F value Pr(>F)
> Residuals 18 381.00   21.17
> --------------------
> 
> What it dont calculate the F and P for treatment?

Would it be easier to do it this way?

> library(lme4)
Loading required package: Matrix
Loading required package: lattice
> (fm1 <- lmer(Glycogen ~ Treatment + (1|Treatment:Rat) + (1|Treatment:Rat:Liver)))
Linear mixed-effects model fit by REML
Formula: Glycogen ~ Treatment + (1 | Treatment:Rat) + (1 | Treatment:Rat:Liver) 
      AIC      BIC    logLik MLdeviance REMLdeviance
 231.6213 241.1224 -109.8106    234.297     219.6213
Random effects:
 Groups              Name        Variance Std.Dev.
 Treatment:Rat:Liver (Intercept) 14.167   3.7639  
 Treatment:Rat       (Intercept) 36.065   6.0054  
 Residual                        21.167   4.6007  
# of obs: 36, groups: Treatment:Rat:Liver, 18; Treatment:Rat, 6

Fixed effects:
            Estimate Std. Error DF t value Pr(>|t|)    
(Intercept) 140.5000     4.7072 33 29.8481   <2e-16
Treatment2   10.5000     6.6569 33  1.5773   0.1243    
Treatment3   -5.3333     6.6569 33 -0.8012   0.4288    
> anova(fm1)
Analysis of Variance Table
          Df  Sum Sq Mean Sq   Denom F value  Pr(>F)  
Treatment  2 123.993  61.996  33.000   2.929 0.06746

The degrees of freedom for the denominator are an upper bound (in this
case a rather gross upper bound) so the p-value is a lower bound.  It
is on my "To Do" list to improve tthis but I have a rather long "To
Do" list.