[R-sig-ME] Anova Random Effect Nested

[Douglas Bates]

On Feb 12, 2008 10:59 AM,  <marcioestat at pop.com.br> wrote:

>  Hi listers,

> I am trying to do an anova with random
> effects and
> my data is nested...
> I didn't find what procedure
> for the random
> effects in anova... The AOV procedure doesn't tell
> about random effects,
> in fact there is no information about nested
> data too... I found the
> function LME, but I am not sure if I can use
> this function... Could
> anybody would send me any example code of an
> anova with random effect and
> nested data...

It is often easier if you could provide the list with a small example
data set and describe the type of model that you want to fit.  The lme
function in the nlme package can fit models with nested random
effects.  Depending on exactly the type of model that you want to fit
the lmer function in the lme4 package can be faster and more reliable
than lme.

One data set that may illustrate fitting nested random effects is the
Oats data in the nlme package (also provided in a somewhat simpler
form in the MEMSS package)

> library(lme4)
Loading required package: Matrix
Loading required package: lattice

Attaching package: 'Matrix'


	The following object(s) are masked from package:stats :

	 xtabs

> data(Oats, package = "MEMSS")
> str(Oats)
'data.frame':	72 obs. of  4 variables:
 $ Block  : Factor w/ 6 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Variety: Factor w/ 3 levels "Golden Rain",..: 3 3 3 3 1 1 1 1 2 2 ...
 $ nitro  : num  0 0.2 0.4 0.6 0 0.2 0.4 0.6 0 0.2 ...
 $ yield  : num  111 130 157 174 117 114 161 141 105 140 ...

The yield is the response for various varieties of oats planted at
different farms (the Block factor) and with different amounts of
nitrogen applied.  In agricultural experimental design we would say
that the farm is a blocking factor, the variety is the whole plot
factor and the level of nitrogen is the subplot factor.  We model
Block and the whole plot as nested random effects.

The enclosed plot of the data, produced by

> xyplot(yield ~ nitro | Block, Oats, groups = Variety, aspect = 'xy',
+        auto.key = list(columns = 3, lines = TRUE), type = "o")

shows whole plot differences but not strong trends for Variety.

To fit a model with Block and Variety within Block as random effects we use

> print(Om1 <- lmer(yield ~ nitro + (1|Block/Variety), Oats))
Linear mixed model fit by REML
Formula: yield ~ nitro + (1 | Block/Variety)
   Data: Oats
 AIC   BIC logLik deviance REMLdev
 603 614.4 -296.5    604.3     593
Random effects:
 Groups        Name        Variance Std.Dev.
 Variety:Block (Intercept) 121.11   11.005
 Block         (Intercept) 210.37   14.504
 Residual                  165.56   12.867
Number of obs: 72, groups: Variety.Block, 18; Block, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept)   81.872      6.945   11.79
nitro         73.667      6.782   10.86

Correlation of Fixed Effects:
      (Intr)
nitro -0.293

The conditional modes (also called the BLUPs) of the random effects are

> ranef(Om1)
An object of class "ranef.mer"
[[1]]
                (Intercept)
Golden Rain:I     3.2331282
Golden Rain:II    4.9719015
Golden Rain:III  -8.0610704
Golden Rain:IV    6.4426777
Golden Rain:V     1.4235656
Golden Rain:VI   -5.6501375
Marvellous:I      0.6246351
Marvellous:II    10.9341715
Marvellous:III   15.6016885
Marvellous:IV    -3.2460109
Marvellous:V     -6.2155928
Marvellous:VI     8.3239327
Victory:I        10.4996447
Victory:II      -14.4054759
Victory:III     -11.2285263
Victory:IV       -5.8545041
Victory:V        -1.1849275
Victory:VI       -6.2091003

[[2]]
    (Intercept)
I     24.939661
II     2.606625
III   -6.406113
IV    -4.616819
V    -10.382321
VI    -6.141032
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Oats.pdf
Type: application/pdf
Size: 21114 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20080212/eb8312c6/attachment.pdf>