How do I make nested x-labels like a JMP variability plot, but using ggplot2

Question

I love the JMP variability plot. (link) It is a powerful tool.

The example the plot has 2 x-axis labels, one for part-number and one for operator.

Here the JMP variability plot displays more than 2 levels of variables. The following splits by oil amount, batch size, and popcorn type. It can take some work to find the right sequence to show strongest separation, but this is an excellent tool for communication of information.

How does one do this, the multiple-level x-labels, with R using the ggplot2 library?

The best that I can find is this (link, link), which separates based on cylinder count, but does not make the x-axis labels.

My example code is this:

#reproducible
set.seed(2372064)

#data (I'm used to reading my own, not using built-in)
data(mtcars)
attach(mtcars)

#impose factors as factors
fact_idx <- c(2,8:11)
for(i in fact_idx){
     mtcars[,i] <- as.factor(mtcars[,i])
}

#boxplot
p <- ggplot(mtcars, aes(gear, mpg, fill=cyl)) + 
     geom_boxplot(notch = TRUE)  
p

The plot this gives is:

How do I make the x-axis lables indicate both gears and cylinders?

In jmp I get this:

Answer 1

You could use R-package VCA which comes with function varPlot implementing variability charts similar to JMP. There are multiple examples provided in the help. Your example would look like this:

library(VCA)
dat <- mtcars[order(mtcars$cyl, mtcars$gear),]

# default
varPlot(mpg~cyl/gear, dat)    
# nicely formatted
varPlot(mpg~cyl/gear, dat, 
        BG=list(var="gear", col=paste0("gray", c(90,80,70)), 
                col.table=T),
        VLine=list(var="cyl"), Mean=NULL,
        MeanLine=list(var=c("cyl", "gear"), col=c("blue", "orange"),
                      lwd=c(2,2)),
        Points=list(pch=16, cex=1))