In pesticide research or, in general, agriculture research, we very commonly encouter experiments with two/three crossed factors and some other treatment that is not included in the factorial structure. For example, let’s consider an experiment with two herbicides (rimsulfuron and dicamba) at two rates (40 and 60 g/ha for rimsulfuron and 0.6 and 1 kg/ha for dicamba) and with four adjuvant treatments (surfactant, frigate, mineral oil and no adjuvant). Apart from this fully crossed structure, we need to introduce, at least, an untreated control and a hand-weeded control. The design for such an experiment has been termed ‘augmented factorial’, because we are, indeed, including some extra treatment levels beyond the crossed factorial structure.

The dataset relating to the above example is named ‘adjuvants.csv’, it is available in an online repository and it can be loaded by using the code below. We have three predictors (Herbicide, Adjuvant and Dose) and two quantitative response variables (WeedCoverage and Yield). We need to transform the variables ‘Herbicide’ and ‘Adjuvant’ into factors, while we recode the two dose levels as ‘HIGH’ and ‘LOW’, in order to make the two herbicides comparable (we could also work with the original dose variable and fit a regression model, but I have made this point in another post, see here).

rm(list = ls())
library(emmeans)
library(multcomp)
fileName <- "https://www.casaonofri.it/_datasets/adjuvants.csv"
dataset <- read.csv(fileName)
head(dataset)
##   Plot   Herbicide   Adjuvant Dose Block WeedCoverage Yield
## 1    1  HandWeeded       None    0     1         0.00 14.55
## 2    2     Dicamba    Frigate    1     1        27.70 11.10
## 3    3     Dicamba Surfactant    1     1        22.60 13.16
## 4    4     Dicamba MineralOil    1     1        16.35 14.22
## 5    5     Dicamba       None    1     1        45.00 14.69
## 6    6 Rimsulfuron MineralOil   60     1        16.65 12.01
dataset$DoseF <- with(dataset, factor(ifelse(Dose==60 | Dose == 1, "HIGH",
                               ifelse(Dose == 0, "NONE", 
                                      "LOW"))))
dataset[,c(1:3,5)] <- lapply(dataset[,c(1:3,5)], function(i) factor(i))

To my experience, such an augmented factorial design may give a few troubles at the data analysis stage. In particular, I have noted that many colleagues, in order to avoid problems, combine the three predictors into one and fit a simple one-way ANOVA model, as shown below, e.g., for the yield response

Comb <- with(dataset, factor(Herbicide:Adjuvant:DoseF))
model1 <- lm(Yield ~ Comb, data = dataset)
anova(model1)
## Analysis of Variance Table
## 
## Response: Yield
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Comb      17 72.604  4.2708  4.8402 2.143e-06 ***
## Residuals 62 54.707  0.8824                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Such an analysis, to my opinion, is rather inefficient. Indeed, the aim of an experiment like this one is not only to compare the treatments with the control, but it is also to analyse the factorial structure, to assess the existence of possible interaction effects among the predictors. Therefore, it may be worth to take a minute to develop a more efficient model to describe our experimental data.

With R, it will be convenient to create a new ‘dummy’ variable to part the controls and the factorial structure, as shown in the box below. In this case, we have two controls, but, in other case, this dummy variable will be binary.

dataset$CvT <- with(dataset, factor(ifelse(Herbicide == "HandWeeded", "chk1",
                      ifelse(Herbicide == "Unweeded", "chk2","trt"))))

Now, we code the model, by nesting the factorial structure within the newly created variable.

model2 <- lm(Yield ~ CvT/(Herbicide*DoseF*Adjuvant), data = dataset)
anova(model2)
## Analysis of Variance Table
## 
## Response: Yield
##                              Df Sum Sq Mean Sq F value    Pr(>F)    
## CvT                           2 46.251 23.1253 26.2080 5.637e-09 ***
## CvT:Herbicide                 1  3.701  3.7008  4.1941  0.044805 *  
## CvT:DoseF                     1  0.106  0.1064  0.1206  0.729531    
## CvT:Adjuvant                  3  6.546  2.1820  2.4729  0.069918 .  
## CvT:Herbicide:DoseF           1  0.220  0.2197  0.2490  0.619534    
## CvT:Herbicide:Adjuvant        3 11.345  3.7815  4.2856  0.008193 ** 
## CvT:DoseF:Adjuvant            3  2.211  0.7372  0.8354  0.479544    
## CvT:Herbicide:DoseF:Adjuvant  3  2.225  0.7415  0.8404  0.476920    
## Residuals                    62 54.707  0.8824                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The ‘CvT’ effects represents the contrast between the two controls and the treated plots and it has three levels (untreated, hand-weeded and treated). We see that this effect is significant, as well as the ‘Herbicide:Adjuvant’ interaction within treated plots. On the other hand, the ‘Dose’ effect within treated plots, is not significant. We can use the emmeans() function to extract the relevant marginal means and compare them.

cld(emmeans(model2, ~CvT), Letters = LETTERS)
##  CvT  emmean    SE df lower.CL upper.CL .group
##  chk2   11.6 0.332 62     11.0     12.3  A    
##  trt    13.9 0.117 62     13.6     14.1   B   
##  chk1   14.8 0.332 62     14.1     15.5    C  
## 
## Results are averaged over the levels of: Adjuvant, DoseF, Herbicide 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
cld(emmeans(model2, ~CvT:Herbicide:Adjuvant), Letters = LETTERS)
##  Herbicide   Adjuvant   CvT  emmean    SE df lower.CL upper.CL .group
##  Unweeded    None       chk2   11.6 0.332 62     11.0     12.3  A    
##  Rimsulfuron None       trt    13.0 0.332 62     12.3     13.6  AB   
##  Dicamba     Frigate    trt    13.3 0.332 62     12.6     13.9   BC  
##  Dicamba     Surfactant trt    13.5 0.332 62     12.8     14.2   BCD 
##  Dicamba     None       trt    13.7 0.332 62     13.1     14.4   BCD 
##  Dicamba     MineralOil trt    14.1 0.332 62     13.4     14.8   BCD 
##  Rimsulfuron MineralOil trt    14.3 0.332 62     13.6     15.0   BCD 
##  Rimsulfuron Surfactant trt    14.6 0.332 62     13.9     15.2    CD 
##  Rimsulfuron Frigate    trt    14.7 0.332 62     14.0     15.3    CD 
##  HandWeeded  None       chk1   14.8 0.332 62     14.1     15.5     D 
## 
## Results are averaged over the levels of: DoseF 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 10 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

I suggest that this type of analysis is simple enough and makes better use of the information contained in our dataset.

Thanks for reading and happy coding!

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