You beat me at RHA, I wrote an almost identical code to you before I saw your message. Anyway, thanks!
I also wanted to remove β2014β and β2015β from the gray boxes (which I did not indicate in my first post), so I had to make some additional changes.
With some inspiration from here , here , and here , I came up with the following (really ugly) code:
data14 <- rnorm(10) data15 <- rnorm(10, mean = 500) A1 <- data.frame("Y"=data14, "X"=1:10, "Q"=1, "year"=2014) A2 <- data.frame("Y"=data14, "X"=1:10, "Q"=2, "year"=2014) A3 <- data.frame("Y"=data14, "X"=1:10, "Q"=3, "year"=2014) A4 <- data.frame("Y"=data14, "X"=1:10, "Q"=4, "year"=2014) N1 <- data.frame("Y"=data15, "X"=1:10, "Q"=1, "year"=2015) N2 <- data.frame("Y"=data15, "X"=1:10, "Q"=2, "year"=2015) N3 <- data.frame("Y"=data15, "X"=1:10, "Q"=3, "year"=2015) N4 <- data.frame("Y"=data15, "X"=1:10, "Q"=4, "year"=2015) A <- rbind(A1, A2, A3, A4) N <- rbind(N1, N2, N3, N4) tmp <- data.frame(rbind(A, N))
Then I made a simple function by naming the variables correctly
labFunc <- function(data, var1, var2, names) { data$id <- NA loop <- length(levels(factor(data[[var2]]))) for (i in 1:loop) { data[data[[var1]] == 2014 & data[[var2]] == levels(factor(data[[var2]]))[i], "id"] <- names[i] data[data[[var1]] == 2015 & data[[var2]] == levels(factor(data[[var2]]))[i], "id"] <- paste(names[i], "") } first <- levels(factor(data$id))[seq(from=1, to = length(levels(factor(data$id))), by = 2)] second <- levels(factor(data$id))[seq(from=2, to = length(levels(factor(data$id))), by = 2)] data$id <- factor(data$id, levels=paste(c(first, second))) return(data) } names <- c("Q1", "Q2", "Q3", "Q4") data <- labFunc(tmp, "year", "Q", names)
Make a plot:
p <-ggplot(data, aes(y = Y, x = X)) + geom_line() + facet_wrap( ~ id , ncol = 4, scales = "free")
And then finally add the main tags
z <- ggplotGrob(p) # New strip at the top z <- gtable_add_rows(z, unit(1, "lines"), pos = 0) # New row added to top z <- gtable_add_rows(z, unit(1, "lines"), pos = 6) # New row added to top #z <- gtable_add_rows(z, unit(9, "lines"), pos = 0) # New row added to top # Check the layout gtable_show_layout(z) # New strip goes into row 2 # New strip spans columns 4 to 8 z <- gtable_add_grob(z, list(rectGrob(gp = gpar(col = NA, fill = grey(0.8), size = .5)), textGrob("2014", vjust = .27, gp = gpar(cex = .75, fontface = 'bold', col = "black"))), 2, 4, 2, 14, name = c("a", "b")) z <- gtable_add_grob(z, list(rectGrob(gp = gpar(col = NA, fill = grey(0.8), size = .5)), textGrob("2015", vjust = .27, gp = gpar(cex = .75, fontface = 'bold', col = "black"))), 7, 4, 7, 14, name = c("a", "b")) # Add small gap between strips - below row 2 z <- gtable_add_rows(z, unit(2/10, "lines"), 2) z <- gtable_add_rows(z, unit(5/10, "lines"), 7) # Draw it grid.newpage() grid.draw(z)

It was a little harder than I thought, but thanks to everyone for their help!