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:

• 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)
• thanks for your help (i have 12h to wait to unlock the +50 :( indeed i was in big doubt concerning the naming of my subplots (4 or 64 unique labels). The figure is correct: irrigation is not "randomized", that's unfortunate i agree (they told me: "too expansive to do it differently"!). Thanks for the links. One more question: i get a residuals plot not looking good: cone-shaped (like this: "<"), error seems proportional to the Y values. is there a way to correct this in this type of model? – agenis Mar 27 '15 at 16:06
• The most obvious solution (and one that often fixes other problems) is to transform the response, most often log-transforming. – Ben Bolker Mar 27 '15 at 16:19