Specifying a model with random effects for a strip-plot designed experiment in R

Question

I am currently trying to understand the analysis of strip-plots (in R) and I came across the example described here (https://www.statforbiology.com/_statbookeng/a-brief-intro-to-mixed-models):

dataset <- read.csv("https://www.casaonofri.it/_datasets/recropS.csv")
head(dataset)
##     Herbicide     Crop Block CropBiomass
## 1       Check soyabean     1    199.0831
## 2       Check soyabean     2    257.3081
## 3       Check soyabean     3    345.5538
## 4       Check soyabean     4    210.8574
## 5 rimsulfuron soyabean     1    225.5651
## 6 rimsulfuron soyabean     2    195.3952
dataset$Herbicide <- factor(dataset$Herbicide)
dataset$Crop <- factor(dataset$Crop)
dataset$Block <- factor(dataset$Block)
dataset$Rows <- factor(dataset$Crop:dataset$Block)
dataset$Columns <- factor(dataset$Herbicide:dataset$Block)

in the strip-plot design, the rows are main plots for columns and vice versa (in analogy to split plots). "each row is uniquely defined by a specific block and crop and each column is uniquely defined by a specific herbicide and block". So far I think I understand.

The model is then fitted as such (still from https://www.statforbiology.com/_statbookeng/a-brief-intro-to-mixed-models):

model.strip <- lmer(CropBiomass ~ Block + Herbicide*Crop + 
    (1|Rows) + (1|Columns), data = dataset)
anova(model.strip, ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Block           21451  7150.3     3 4.1367  2.5076 0.19387  
## Herbicide         148   147.9     1 3.0000  0.0519 0.83450  
## Crop            43874 21936.9     2 6.0000  7.6932 0.02208 *
## Herbicide:Crop  12549  6274.4     2 6.0000  2.2004 0.19198  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1. In the model specification ~Block + Herbicide * Crop + (1|Rows) + (1|Columns), data = dataset) the block is stated as a fixed effect. As I percieve it, the block is already considered for in the Row and in the Column variable. Is it only specified in the fixed effects, because there is also a special interest in the specific differences between blocks? And if I would consider the block as a raondom effect, would it just add the term + (1|Block) to the model? i.e. `lmer(CropBiomass ~ Herbicide*Crop + (1|Rows) + (1|Columns) + (1|Block), data = dataset)

2. If there would be an additional fixed effect (e.g. repeating the experiment exactly the same way for 3 years, considering the years as fixed (factorial) effect), would the model then be CropBiomass ~ Herbicide * Crop * Year + (1|Rows) + (1|Columns) + (1|Block)

3. If there were subsamples/pseudoreplicates in each plot (e.g. a plot with specific Herbicide application and specific Crop has 2 plants, and the biomass of each individual plant is measured), how would the model look like? would it be CropBiomass ~ Herbicide * Crop * Year + (1|Rows) + (1|Columns) + (1|Block:PlotID)

I am glad for any help to further understand the analysis of split plots, thank you!

Answer 1

Let's name the models under consideration for ease of reference. Using R's formula notation:

M1. CropBiomass ~ Herbicide * Crop + Block + (1 | Row) + (1 | Column)
M2. CropBiomass ~ Herbicide * Crop + (1 | Block) + (1 | Row) + (1 | Column)

It's easy to check that on the example data M1 and M2 give the same estimates for the fixed effects of Herbicide, Crop & their interaction, which are presumably the focus of the analysis. The reason is that the experiment is balanced. So — as long as there is balance (no data point is lost) — it doesn't matter whether the Block effect is random or fixed.

In the split-plot design rows and columns are nested within blocks. So I would use instead:

M3. CropBiomass ~ Herbicide * Crop + (1 | Block) + (1 | Block:Row) + (1 | Block:Column)

Again, it's easy to check that on the example data M2 and M3 give exactly the same estimates for both random and fixed effects. The reason is that the rows and columns have unique IDs across blocks.

For example, here are the plots in the first row:

#>   Block Row        Column
#>   1     1.soyabean 1.Check      
#>   2     2.soyabean 2.Check      
#>   3     3.soyabean 3.Check      
#>   4     4.soyabean 4.Check      
#>   1     1.soyabean 1.rimsulfuron
#>   2     2.soyabean 2.rimsulfuron
#>   3     3.soyabean 3.rimsulfuron
#>   4     4.soyabean 4.rimsulfuron

If instead rows were labeled 1,2,3 and columns 1,2 in each of the blocks (so that the same labels 1,2,3 appear across blocks), then M2 and M3 wouldn't specify the same model. Try it out by re-indexing rows and columns such that the plots in the first row are identified by:

#>   Block Row   Column
#>   1     1     1      
#>   2     1     1      
#>   3     1     1      
#>   4     1     1      
#>   1     1     2      
#>   2     1     2      
#>   3     1     2      
#>   4     1     2

It's preferable that models to not depend on the vagaries of ID labeling, so I'd use M3.

To add a time component and analyze experimental data collected over 3 years, I would start with the model:

CropBiomass ~ factor(Year) * Herbicide * Crop + (1 | Block) + (1 | Block:Row) + (1 | Block:Column)

Or perhaps the model with all two-way product terms but without the three-way product term:

CropBiomass ~ (factor(Year) + Herbicide + Crop)^2 + (1 | Block) + (1 | Block:Row) + (1 | Block:Column)

I assume that the multi-year experiment will be performed on the same blocks, with the assignment of crops & herbicides to plots within blocks randomized each year. In that case it's likely that 1.soyabean in Year 1 won't indicate the same plot as 1.soyabean in Years 2 and 3. So for the multi-year experiment it would be better to assign row and column IDs explicitly rather than do it through a shortcut such as dataset$Crop:dataset$Block.

Finally, if there are two or more plants measured per plot, as suggested, a random plot-within-block effect will allow for correlation between the biomass of plants in the same plot.

M4. CropBiomass ~ Herbicide * Crop + (1 | Block) + (1 | Block:Row) + (1 | Block:Column) + (1 | Block:Plot)

If the design is balanced so that there are k plants in each plot/block/year, it's also valid to average the biomass of the k plants and use M3. This analysis will give the same inference for the herbicide and crop effects. (The experimental units are the plots, not the plants. Averaging the biomass of plants in a plot means that the biomass of an experimental unit is observed with less error. As long as we average the same number of plants, the error is reduced by the same factor in every plot.)

PS: Don't use Type III sums of squares. Here is why: How to interpret type I, type II, and type III ANOVA and MANOVA?.