[R] How to run Linear mixed model for an experiment design with species nested in an random block experiment

[Faming Wang]

Dear all,

   I  have conducted an N and P field addition experiment in a tropical
forest, and we used a random block design in this experiment, briefly, we
had four plots in each block (Control, +N， +P, and +NP), and five blocks
located in the forest randomly. Totally we have 20 plots, with four
treatments and five replicated blocks. In each plot, we selected five
species  plants (some plots only contains 3 or 4 species) to measure their
leave variables, like N concentration.  We want to know the effect of N and
P addition as well as the species level variability (inter-species )  of
leaf N.  So we used linear mixed effect models to conduct our statistical
analysis, the sample code was listed below. Can anybody take a look at this
script, and help me to figure out how to analysis species effect using LME?

 Thanks!


R script attached

####leaf N concentration

### FIRST, WE TEST WHETHER NESTING SPECIES AS A RANDOM EFFECT IMPROVES
THE FULL MODEL

lmeleafN1<-lme(fixed=N~Naddition*Paddition*Species, random = ~1|Block,
data = NPdata, method = "ML", na.action=na.exclude)
lmeleafN1a<-lme(fixed=N~Naddition*Paddition*Species, random =
~1|Block/Species, data = NPdata, method = "ML", na.action=na.exclude)
anova(lmeleafN1, lmeleafN1a )

### NESTING SPECIES WITHIN BLOCK DOESN'T IMPROVE THE MODEL, SO WE CAN
USE THE MODEL WITH THE SIMPLER RANDOM EFFECT

lmeleafN2<-lme(fixed=N~Naddition*Paddition+Species, random = ~1|Block,
data = NPdata, method = "ML", na.action=na.exclude)
lmeleafN3<-lme(fixed=N~Naddition+Paddition*Species, random = ~1|Block,
data =NPdata, method = "ML", na.action=na.exclude)
lmeleafN4<-lme(fixed=N~Naddition+Paddition+Species, random = ~1|Block,
data = NPdata, method = "ML", na.action=na.exclude)

AIC(lmeleafN1, lmeleafN2, lmeleafN3, lmeleafN4)

# THE FULL MODEL CLEARLY HAS THE LOWEST AIC
# CHECK AGAINST THE NULL MODEL
lmeleafN0<-lme(fixed=N~1, random = ~1|Block, data = NPdata, method =
"ML", na.action=na.exclude)

anova(lmeleafN1, lmeleafN0)

## AND CHECK THE MODEL FIT WITH DIAGNOSTIC PLOTS

par(mfrow = c(2,2))
plot(resid(lmeleafN1)  ~ fitted(lmeleafN1))
abline(h=0, lty=2)
hist(resid(lmeleafN1))
qqnorm(resid(lmeleafN1))
qqline(resid(lmeleafN1))
anova(lmeleafN1)



-- 




Sincerely

Faming Wang

Associate Scientist
Deputy Director of Xiaoliang Research Station,
South China Botanical Garden, Chinese Academy of Sciences
Xingke Road 723, Guangzhou, China. 519650
Email: wangfm at scbg.ac.cn
Tel/Fax:0086-20-37252905

	[[alternative HTML version deleted]]