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

# Reference

- Piepho, H.P., Edmondson, R.N., 2018. A tutorial on the statistical analysis of factorial experiments with qualitative and quantitative treatment factor levels. J Agronomy Crop Science 204, 429–455. https://doi.org/10.1111/jac.12267