A question on using LME model in species nested within a random block design

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 randomly distributed plots in each block (Control, +N? +P,
and +NP), and five blocks located in the forest. Totally we have 20 plots,
with two N treatments and two P treatment 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, P
concentration, and photosynthesis rate et al.  We want to know the effect
of N and P addition as well as the species level changes (inter-species )
 on leaf variables. Since some plots some specific species are missing in
some plots some specific species, it was unbalanced at the species level.
We used linear mixed effect models to conduct our statistical analysis:

  We firstly tested the random effect with blocks, and species within plots
within blocks, and found that nesting plots and species within block did
not improve the model fitness, so we choose only block as random effect.
For fixed effects, N-addition, P-addition, species and their interaction
were considered fixed effects in models. The significance of each term was
determined by comparing nested models using likelihood ratio tests and AICs
to check for model improvement. Since there was better model fit (lower AIC
values) with interaction terms, we selected the full factor model. However,
as there was a highly significant effect of tree species identity and
species related interactions, species-specific responses to N- and
P-addition were also investigated with separate models with N, P and their
interaction as fixed effects and block as a random effect.
 however, our reviewers were not happy with this statistical methods and
pointed out that "Species is treated as a fixed factor, generating a
three-way factorial ANOVA. Species cannot be treated in this way because
all five species were present in each 10-by-10-m plot. To implement a
three-way ANOVA design, the entire experiment (five blocks of the four
factorial N and P treatments) would have to be repeated once for each
species. Species cannot be treated as a fixed factor because all five
species were measured in the same experimental plots. This is a split-plot
design. Alternatively, MANOVA might be performed treating the five species
as five response variables. A split-plot design or a MANOVA approach would
allow the authors to investigate interspecific variation in responses."

  I am very confused on the reviewer's comments,  it seems to me that the
reviewer compared our LME model with 3-way ANOVA. If we used 3-way ANOVA, I
know that my experiment is species nested in a random block design, and we
could not directly use 3-Way ANOVA, which the error df would be
overestimated.

   Below I attached my sample data and my current R script for LME model in
dropbox. See below links:
https://www.dropbox.com/s/6kd3kq5mlyuyqz6/NPrawdata.csv?dl=0

https://www.dropbox.com/s/fpqdbm6go0g8ak0/Faming%20NP%20model.R?dl=0


--

Sincerely

Faming Wang

        [[alternative HTML version deleted]]