[R] Trying to make a nested lme analysis

[Douglas Bates]

Where is the rats data available?

It looks as if you have an lme model with both a fixed effect for
Treatment and a random effect for Treatment.  I would guess that you
want to have a fixed effect for treatment and random effects for

 Rat %in% Treatment 

and 

 Liver %in% Rat %in% Treatment

If so you would first create a factor for Rat %in% Treatment, say rTrT
by 

 rats$rTrt = getGroups(~ 1 | Treatment/Rat, data = rats, level = 2)

then fit the lme model as

 lme(Glycogen ~ Treatment, data = rats, random = ~ 1|rTrT/Liver)


"Ronaldo Reis Jr." <chrysopa at insecta.ufv.br> writes:

> Hi,
> 
> I'm trying to understand the lme output and procedure.
> I'm using the Crawley's book.
> 
> I'm try to analyse the rats example take from Sokal and Rohlf (1995).
> I make a nested analysis using aov following the book.
> 
> > summary(rats)
>     Glycogen       Treatment      Rat          Liver  
>  Min.   :125.0   Min.   :1   Min.   :1.0   Min.   :1  
>  1st Qu.:135.8   1st Qu.:1   1st Qu.:1.0   1st Qu.:1  
>  Median :141.0   Median :2   Median :1.5   Median :2  
>  Mean   :142.2   Mean   :2   Mean   :1.5   Mean   :2  
>  3rd Qu.:150.0   3rd Qu.:3   3rd Qu.:2.0   3rd Qu.:3  
>  Max.   :162.0   Max.   :3   Max.   :2.0   Max.   :3  
>
> > attach(rats)
> > Treatment <- factor(Treatment)
> > Rat <- factor(Rat)
> > Liver <- factor(Liver)
> 
> > 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               
> > 
> 
> OK,
> 
> Then I try to make this analysis using lme.
> 
> > model <- lme(Glycogen~Treatment, random=~1|Treatment/Rat/Liver)
> > summary(model)
> Linear mixed-effects model fit by REML
>  Data: NULL 
>        AIC      BIC    logLik
>   233.6213 244.0968 -109.8106
> 
> Random effects:
>  Formula: ~1 | Treatment
>         (Intercept)
> StdDev:    3.541272
> 
>  Formula: ~1 | Rat %in% Treatment
>         (Intercept)
> StdDev:     6.00658
> 
>  Formula: ~1 | Liver %in% Rat %in% Treatment
>         (Intercept) Residual
> StdDev:    3.764883 4.600247
> 
> Fixed effects: Glycogen ~ Treatment 
> Error in if (any(wchLv <- (as.double(levels(xtTab[, wchPval])) == 0))) { : 
> 	missing value where logical needed
> In addition: Warning message: 
> NaNs produced in: pt(q, df, lower.tail, log.p) 
> > 
> 
> The random effects are correct, the variance component is OK:
> 
> In nested aov | In nested lme
> Residual
> 21.1666       | 21.16227
> Liver in Rats
> 14.16667      | 14.17434
> Rats in Treatment
> 36.0648       | 36.079
> 
> But I not understand why the Fixed effects error?
> 
> What is the problem in my formula to make this analysis using lme?
> 
> Thanks for all
> Inte
> Ronaldo
> -- 
> Anger kills as surely as the other vices.
> --
> |   // | \\   [*****************************][*******************]
> || ( õ   õ )  [Ronaldo Reis Júnior          ][PentiumIII-600     ]
> |      V      [UFV/DBA-Entomologia          ][HD: 30 + 10 Gb     ]
> ||  /     \   [36571-000 Viçosa - MG        ][RAM: 128 Mb        ]
> |  /(.''`.)\  [Fone: 31-3899-2532           ][Video: SiS620-8Mb  ]
> ||/(: :'  :)\ [chrysopa at insecta.ufv.br      ][Modem: Pctel-onboar]
> |/ (`. `'` ) \[ICQ#: 5692561                ][Kernel: 2.4.18     ]
> ||  ( `-  )   [*****************************][*******************]
> ||| _/   \_Powered by GNU/Debian W/Sarge D+ || Lxuser#: 205366
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help

-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/