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:

[hidden email]
Tel/Fax:0086-20-37252905

[[alternative HTML version deleted]]

______________________________________________

[hidden email] mailing list -- To UNSUBSCRIBE and more, see

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.