[R] aov vs lme for split plot analysis

[Rafael Rubio de Casas]

Dear R community,

I am trying to do a split-plot analysis as follows. I have a data set 
(“morf”) with plant data from 6 “blocks” at different latitudes, each 
divided in 3 plots. The full-plot “treatment” is “soil type” and has 
three levels. Within each plot I have two levels of radiation, coded as 
“SUN” and “SHADE”. I have data for several response traits for 30 plants 
within each subplot, that is 1080 plants in total (6 Blocks x 3Plots x 
2Radiation levels x 30 Plants).
I want to measure a) if there are differences between latitudes 
(“blocks”) and b) how these differences interact with radiation levels. 
I am not so interested in the effect of soil type, which I already know, 
but radiation effects depend on soil type. The model I’m looking for 
should look something like:
Response ~ Block*Radiation | Soil type
I have tried to use lme. The basic formulation I have tried is:
lme(Resp~Block*Radiation, random=~ 1 | Plot/Radiation, data=morf)
but I keep on getting the error message:
“Error in getGroups.data.frame(dataMix, groups) : Invalid formula for 
groups”
I don’t know why this is happening, could it be because “Radiation” has 
only two levels? I don't have missing values or any other obvious source 
of "noise"
Since neither lme nor lmer weren’t working, I tried “aov”. I played 
around a bit with the "Oats" data set and saw that:
aov( yield ~ ordered(nitro) * Variety + Error(Block/Variety), data = Oats)
gives the same results as the standard
lme( yield ~ ordered(nitro) * Variety, data = Oats, random = ~ 1 | 
Block/Variety )
So I used the following formulation
aov(Resp~Block*Radiation+Error(Plot/Radiation), data=morf)
This seems to work fine, but I am not that confident I am using the 
correct syntax.
Any advice and/or guidance would be much appreciated.

Thank you in advance.

RAFA