# How to code a two-way ANOVA with split-plot design using R?

I am analyzing a split-plot land-use experiment and I am stuck on the R code.

We used four sites and on each site the following design: We conducted a full-factorial design with two factors (land-use intensity and land-use method) with two levels, resulting in four plots per site. Then we divided each plot in three and performed a soil treatment with three levels (resulting in 12 subplots per site).

This is my first split-plot analysis and maybe I am just confused, because I read through too many examples and varying R codes. However, I can´t find a reference or code suitable for my experiment, meaning a split-plot analysis with a two-way anova on top. For example, I am having difficulties to assign an order in the error term for the land-use factors.

This is my first trial so far (heavily based on examples by Dormann 2017) and I am not happy with it:

lme(x~intensity*method*soil, random=~1|site/soil)


I am grateful for any suggestions and corrections.