Motivating example

Let’s think of a field experiment, where 6 genotypes of faba bean were compared under two different sowing times (autumn and spring). For the ease of organisation, such an experiment was laid down as a split-plot, in a randomised complete block design; sowing times were randomly allocated to main-plots, while genotypes were randomly allocated to sub-plots. Let’s stop for a moment… does this sound strange to you? Do you need further information about split-plot designs? You can get some general information at this link and hints on how to analyse the results at this other link

The above experiment was repeated in three years and two locations (six environments in all), in order to explore the environmental variability of results (we will not make any distinction between years and locations, for the sake of this post). In the end, we recorded crop yield and produced a dataset with 288 record (6 environments by 2 sowing times by 6 genotypes by 4 blocks). If you are interested in more detail about this experiment, you can find them in Stagnari et al., (2007).

The resulting dataset (‘fabaBeanSplitMet.csv’) is available in a public online repository and contains six columns, the ‘Location’, the ‘Year’, the ‘Sowing’ time, the ‘Genotype’, the ‘Block’ and the response variable, i.e. the ‘Yield’. After loading the dataset, we need to recode the independent variables into factors and create the new ‘Environment’ factor, as the combination of ‘Year’ and ‘Location’ levels. In the box below, we use the ‘dplyr’ package to accomplish this preliminary step (Wickham et al., 2022).

library(tidyverse)
rm(list=ls())
fileName <- "https://www.casaonofri.it/_datasets/fabaBeanSplitMet.csv"
dataset <- read.csv(fileName)
dataset <- dataset %>% 
  mutate(across(c(Location, Year, Sowing, Genotype, Block),
                .fns = factor),
         Environment = factor(Location:Year))
head(dataset)
##   Location      Year Sowing  Genotype Block Yield       Environment
## 1  papiano 2002-2003 autumn    Chiaro     1  2.05 papiano:2002-2003
## 2  papiano 2002-2003 autumn    Chiaro     2  2.50 papiano:2002-2003
## 3  papiano 2002-2003 autumn    Chiaro     3  2.64 papiano:2002-2003
## 4  papiano 2002-2003 autumn    Chiaro     4  2.45 papiano:2002-2003
## 5  papiano 2002-2003 autumn Collameno     1  2.01 papiano:2002-2003
## 6  papiano 2002-2003 autumn Collameno     2  2.19 papiano:2002-2003

Building a valid model

A model is identified by listing all the effects which we need to explain the observed yield. In this case, considering the aims of our experiment, it is pretty easy to grasp the importance of the ‘sowing date’ effect, the ‘genotype’ effect and their interaction. These are the so-called treatment factors and we have no doubt that they should be included in our model. Furthermore, we should also be interested to know whether those treatment effects interact with the environment effect, so we should clearly add to the model the ‘sowing time by environment’, ‘genotype by environment’ and ‘sowing time by genotype by environment’ interactions.

At this step, it is possible that we have no specific interest in any other effects, apart from those we have just mentioned; however, if we stop now, our model is still incomplete and, therefore, invalid. Indeed, we should also think about possible grouping factors. You may wonder: what are the grouping factors? This aspect needs particular attention.

In split-plot and other very common types of designs, the experimental units are not completely randomised, but they are organised (‘grouped’, indeed) by way of some innate attribute, such as the environment or block or main-plot, which they belong to. These attributes are known as ‘grouping factors (see Piepho et al., 2003) and they introduce a sort of ’parenthood’, so that some observations are more alike than others, because they belong to the same ‘group’ (e.g., same block or same main-plot). If we neglect the effects of ‘grouping factors’, these ‘parenthood’ effects remain in the residuals, which will be correlated. The correlation of residuals represent an important violation of the basic assumptions for linear model fitting and, therefore, the model will be invalid and our paper will be rejected. One first conclusion: please, do never forget the grouping factors, if you want your paper to be accepted!.

What are the grouping factors in this case? First of all we have the environments (six levels), then we have the blocks within each environment (24 levels in all) and, finally, we have the main-plots within each block and within each environment (48 levels, in all). In this latter respect, we can see that each main-plot can be uniquely identified by the combination of one environment, one block and one sowing time. Consequently, the final (valid) model is:

Yield ~ Environment + Sowing + Environment:Sowing + Genotype + 
        Environment:Genotype + Environment:Sowing:Genotype + 
        Environment:Block + Environment:Block:Sowing

In R, we can abbreviate as:

Yield ~ Environment * Sowing * Genotype + 
        Environment:Block + Environment:Block:Sowing

Sorry, I know I am running the risk of being regarded as a boring professor; but, please, remember: failing to include any of the above mentioned effects in the model, unless they are clearly non-significant, leads to totally invalid results!

Now, we need to take a very important decision: which factors are fixed and which factors are random? The rule is that all factors that reference randomisation units (to which treatments are allocated) NEED TO BE RANDOM, while, for the other factors, we can make our own subjective choice. Here, the main-plot factor, to which we allocated the sowing dates, needs to be taken as random. For the other factors, we make the most traditional choice of taking them as fixed, although we need to consider that, in other instances, it might be appropriate to regard the ‘environment’ and ‘block’ effects as random (relating to block effects, you may read Dixon, 2016, for interesting information).

If I were to suggest a simple package to fit the above model, I’d say that you should favour the lmer() function in the lme4 package, where the random effects are coded by using the ‘(1|effect)’ notation, as shown in the box below; before fitting, we load the ‘lme4’ package, together with the ‘lmerTest’ package, which gives us extra-flexibility to produce an ANOVA table:

library(lme4)
library(lmerTest)
modMix <- lmer(Yield ~ Environment * Sowing * Genotype +
               Environment:Block + (1|Environment:Block:Sowing),
               data = dataset)
anova(modMix)
## Type III Analysis of Variance Table with Satterthwaite's method
##                             Sum Sq Mean Sq NumDF DenDF  F value    Pr(>F)    
## Environment                 71.084  14.217     5    18 125.4463 2.430e-13 ***
## Sowing                      45.947  45.947     1    18 405.4282 8.574e-14 ***
## Genotype                     8.030   1.606     5   180  14.1707 1.091e-11 ***
## Environment:Sowing          11.022   2.204     5    18  19.4520 1.086e-06 ***
## Environment:Genotype         9.468   0.379    25   180   3.3418 1.454e-06 ***
## Sowing:Genotype              5.340   1.068     5   180   9.4231 5.388e-08 ***
## Environment:Block            4.398   0.244    18    18   2.1560    0.0561 .  
## Environment:Sowing:Genotype  7.912   0.316    25   180   2.7925 4.513e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Obviously, multiple comparison testing can be done with the ‘emmeans’ package as we have shown elsewhere. Transforming the environment or block effects into random effects is pretty straightforward, by changing the R notation; please remember that, if you regard the environment as random, all its interactions should also be regarded most often regardes as well as random.

Did I menage to make myself clear? If not, drop me a line to the address below. Happy coding!

Prof. Andrea Onofri
Department of Agricultural, Food and Environmental Sciences
University of Perugia (Italy)
Send comments to: andrea.onofri@unipg.it

References

1. Dixon, P., 2016. Should blocks be fixed or random? Conference on Applied Statistics in Agriculture. https://doi.org/10.4148/2475-7772.1474
2. Piepho, H.-P., Büchse, A., Emrich, K., 2003. A Hitchhiker’s Guide to Mixed Models for Randomized Experiments. Journal of Agronomy and Crop Science 189, 310–322.
3. Stagnari, F., Onofri, A., Jemison, J.J., Monotti, M., 2007. Improved multivariate analyses to discriminate the behaviour of faba bean varieties. Agronomy For Sustainable Development 27, 387–397.
4. Wickham H, François R, Henry L, Müller K (2022). Dplyr: A Grammar of Data Manipulation. R package version 1.0.9, https://CRAN.R-project.org/package=dplyr.