[R] how to perform 5 level nested anova

Dieter Menne dieter.menne at menne-biomed.de
Wed Feb 11 08:41:41 CET 2009


Mao Jianfeng <jianfeng.mao <at> gmail.com> writes:

> I am new to R. And, I want to perform a multiple nested anova on a large
> datasets (with 9448 observations).
> 
> I tried the function ("lmer()" in package "lme4"). But, I failed. Can anyone
> help me?
> 
> my dataset("SeedL.txt") was attached. Data are not sorted by factors.
> 
> In this dataset, "SpecN" "PopN"  "TreeN" "ConeN" "SeedN" were 5 factors (as
> Explanatory), with "PopN" nested within "SpecN"; "TreeN" nested within
> "PopN"; "ConeN" nested within "TreeN" and "SeedN" nested within "ConeN".
> 
> "SeedL" is a dependent variate (as Response).
> 
> The model I have ever used in R is:
> 
> cmod<-lmer(SeedL~1+(1|SpecN)+(1|SpecN:PopN)+(1|SpecN:PopN:TreeN),
> data=seedL)
> 
> What I got in R is:
> 
> error: length(f1) == length(f2) is not TRUE
> In addition: Warning messages:
> 1: In SpecN:PopN : 数值表达式一共有9447元素: 只用了第一个 (In Enlish: there
are 9447
> elements in the expression, but only the first one has been used.)

First, note that one of the data was missing in line 51xx, this could have 
cause you some problem.

You probably forgot to convert the integer variables to factors,
these must be marked as categorical variables.
Your model works now, but I doubt you really wanted a random-only
model. One of my guesses, assuming SpecN is the fixed variable,
is shown below. 

Also check if package nlme; it should be possible to handle the
model with that package, which is much better documented
(Book by Pinheiro/Bates).

Dieter



seedL = read.table("seed.txt",header=TRUE)
# Make factors
for (i in 1:5) seedL[,i] = as.factor(seedL[,1])

library(lme4)
# works now, but I doubt you really want a random-only model
cmod<-lmer(SeedL~1+(1|SpecN)+(1|SpecN:PopN)+(1|SpecN:PopN:TreeN),
           data=seedL)
summary(cmod)
# My guess: this is closer to what you want, testing species
cmod1<-lmer(SeedL~SpecN+(1|SpecN:PopN)+(1|SpecN:PopN:TreeN),
           data=seedL)
summary(cmod1)




More information about the R-help mailing list