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.

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:

• 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)