Mixed-effect model design with a sampling variable I am trying to specify a formula for a linear mixed effect model (with lme4) for my experimental design, but I'm not sure I'm doing it right.
The design: basically I'm measuring a response parameter on plants. I have 4 levels of treatment, and 2 irrigation levels.  The plants are grouped in 16 plots, within each plot I sample 4 sub-plots. In each sub-plot I take between 15 and 30 observations (depending on the number of plants found). That is, there are a total of 1500 rows. 

Initially the subplot level was just here for sampling purposes, but I thought I'd like to take it into account in the model (as a 64-level variable) because I saw there was a lot of variability from one sub-plot to another, even inside the same plot (greater than the variability between whole plots).
My first idea was to write:
library(lme4)
fit <- lmer(y ~ treatment*irrigation + (1|subplot/plot), data=mydata)

or 
fit <- lmer(y ~ treatment*irrigation + (1|subplot) + (1|plot), data=mydata)

Is that correct? I'm not sure if I must keep both plot/subplot levels in my formula. 
No fixed effect is significant but the random effects are very significant.
 A: Your model should be written as
fit <- lmer(y ~ treatment*irrigation + (1|plot/subplot), data=mydata)

as subplots are nested within site. although (1|plot) + (1|subplot) would work if the subplots are uniquely labeled (i.e. 1A,1B,1C,...,2A,2B,2C rather than A,B,C...,A,B,C).  My book chapter from Fox et al. Ecological Statistics describes an example of nesting:

On the other hand, in the tick example each chick occurs in only one brood,
  and each brood occurs in only one site: the model specification is (1 | SITE/BROOD/INDEX), read as “chick (INDEX) nested within brood nested within site,” or equivalently (1 | SITE) + (1 | SITE:BROOD) + (1 | SITE:BROOD:INDEX). If the broods and chicks are uniquely labeled, so that the software can detect the nesting, (1 | SITE) + (1 | BROOD) + (1 | INDEX) will also work (do not use (1 | SITE) + (1 | SITE/BROOD) + (1 | SITE/BROOD/INDEX); it will lead to redundant terms in the model).

Other thoughts:


*

*more information on nesting and model specifications at http://glmm.wikidot.com/faq

*are your irrigation treatments really organized as shown in the schematic above, i.e. non-interspersed? Or is that just for convenience of graphical presentation?  If the former, then you have a potentially problematic experimental design ...

*Since subplots are nested within sites, it would be just fine inferentially (following Murtaugh 2007 Ecology "Simplicity and complexity in ecological data analysis") to take the plot means and analyze the data at the plot level.

*For what it's worth, I think you could go even farther and aggregate to the plot level; then you could skip mixed models entirely and just do lm(y~treatment*irrigation, data=my_aggregated_data)
