It's possible to easily do with base graphics :
data <- data.frame(x = rep(LETTERS[1:2], each = 10), y = rep(LETTERS[3:4], 10), z = rnorm(20))
boxplot(z ~ x:y, data = data, axes = F, col = rep(c("red", "blue"), each = 2))
mtext(c("A", "B"), 1, 0, at = c(1.5, 3.5))
mtext(rep(c("C", "D"), 2), 1, 1, at = 1:4)
mtext(expression(bold(" Region")), 1, 0, at = 0.2)
mtext(expression(bold("Treatment")), 1, 1, at = 0.2)
axis(1, col = "black", labels = F)
axis(2)
How can I do the same with ggplot2 boxplots ? I don't want the factor y to be a boxplot colour, as previous solutions have done. In fact, I intend the boxes to be coloured by factor x.