I am trying to combine facet strips across two adjacent panels (there is always two adjacent ones with the same first ID variable, but with two different scenarios, let's call them "A" and "B"). I am not particularly wedded to the gtable + grid solution I tried, but sadly I cannot use the facet_nested() from the ggh4x package (I cannot install it on my company's server due to various restrictions that are in place and needed dependencies - I looked at using only the relevant code, but that again is not easy due to the dependencies).

A minimum viable example of the basic plot I want to make easier to read by indicating which panels "belong together" by combining the top facet strips looks like this:

library(tidyverse)
library(gtable)
library(grid)

idx = 1:16

p1 = expand_grid(id=idx, id2=c("A", "B"), x=1:10) %>%
  mutate(y=rnorm(n=n())) %>%
  ggplot(aes(x=x,y=y)) +
  geom_jitter() +
  facet_wrap(~id + id2, nrow = 4, ncol=8)

The strips with the "1"s, the ones with the "2"s etc. should be combined (in reality it's a somewhat longer text, but this is just for illustration). I was trying to adapt an answer for a similar scenario (https://stackoverflow.com/a/40316170/7744356 - thank you @markus for finding it again), but this is what I tried. As you can see below, the height of what I produce seems wrong. I assume this must be some trivial thing I am overlooking/not understanding.

# Combine strips for a ID
g <- ggplot_gtable(ggplot_build(p1))
strip <- gtable_filter(g, "strip-t", trim = FALSE)
stript <- which(grepl('strip-t', g$layout$name))
  
stript2 = stript[idx*2-1]
top <- strip$layout$t[idx*2-1]
# # Using the $b below instead of b = top[i]+1, also seems  not to work
#bot <- strip$layout$b[idx*2-1] 
l   <- strip$layout$l[idx*2-1]
r   <- strip$layout$r[idx*2]
  
mat   <- matrix(vector("list",
                       length = length(idx)*3),
                nrow = length(idx))
mat[] <- list(zeroGrob())

res <- gtable_matrix("toprow", mat,
                     unit(c(1, 0, 1), "null"),
                     unit( rep(1, length(idx)),
                           "null"))

for (i in 1:length(stript2)){
  if (i==1){
    zz <- res %>% 
      gtable_add_grob(g$grobs[[stript2[i]]]$grobs[[1]], 1, 1, 1, 3) %>%
      gtable_add_grob(g, ., 
                      t = top[i],  
                      l = l[i],  
                      b = top[i]+1,  
                      r = r[i], 
                      name = c("add-strip")) 
  } else {
    zz <- res %>% 
      gtable_add_grob(g$grobs[[stript2[i]]]$grobs[[1]], 1, 1, 1, 3) %>%
      gtable_add_grob(zz, ., 
                      t = top[i],  
                      l = l[i],  
                      b = top[i]+1,  
                      r = r[i], 
                      name = c("add-strip"))
  } 
}

grid::grid.draw(zz)

------------ Update with a ggh4x implementation -----------------

This may solve this type of problem for many, but has its downsides (e.g. axes alignment across rows gets a bit manual, probably need to manually remove x-axes and ensure the limits are the same, add a unified y-axis label, requires installation of a package from github: devtools::install_github("teunbrand/ggh4x@v0.1") for a specific version, plus cowplot interacts badly with e.g. ggtern). So I'd love it, if someone still managed to do a pure gtable + grid version.

library(tidyverse)
library(ggh4x)
library(cowplot)

plots = expand_grid(id=idx, id2=c("A", "B"), x=1:10) %>%
  mutate(y=rnorm(n=n()),
         plotrow=(id-1)%/%4+1) %>%
  group_by(plotrow) %>%
  group_map( ~ ggplot(data=.,
                      aes(x=x,y=y)) +
               geom_jitter() +
               facet_nested( ~ id + id2, ))
            
plot_grid(plotlist = plots, nrow = 4, ncol=1)

Here's a reprex of a somewhat pedestrian way to do it in grid. I have made the "parent" facet somewhat darker to emphasise the nesting, but if you prefer the color to match just change the rectGrob fill color to "gray85".

# Set up plot as per example

library(tidyverse)
library(gtable)
library(grid)

idx = 1:16

p1 = expand_grid(id=idx, id2=c("A", "B"), x=1:10) %>%
  mutate(y=rnorm(n=n())) %>%
  ggplot(aes(x=x,y=y)) +
  geom_jitter() +
  facet_wrap(~id + id2, nrow = 4, ncol=8)

g <- ggplot_gtable(ggplot_build(p1))

# Code to produce facet strips

stript <- grep("strip", g$layout$name)

grid_cols <- sort(unique(g$layout[stript,]$l))
t_vals <- rep(sort(unique(g$layout[stript,]$t)), each = length(grid_cols)/2)
l_vals <- rep(grid_cols[seq_along(grid_cols) %% 2 == 1], length = length(t_vals))
r_vals <- rep(grid_cols[seq_along(grid_cols) %% 2 == 0], length = length(t_vals))
labs   <- levels(as.factor(p1$data$id))

for(i in seq_along(labs))
{
  filler <- rectGrob(y = 0.7, height = 0.6, gp = gpar(fill = "gray80", col = NA))
  tg    <- textGrob(label = labs[i], y = 0.75, gp = gpar(cex = 0.8))
  g     <- gtable_add_grob(g, filler, t = t_vals[i], l = l_vals[i], r = r_vals[i], 
                           name = paste0("filler", i))
  g     <- gtable_add_grob(g, tg, t = t_vals[i], l = l_vals[i], r = r_vals[i], 
                           name = paste0("textlab", i))
}

grid.newpage()
grid.draw(g)

And to demonstrate changing the rectGrob to 50% height and "gray85":

Or if you wanted you could assign a different fill for each cycle of the loop:

Obviously the above method might take a few tweaks to fit other plots with different numbers of levels etc.

^{Created on 2020-07-04 by the reprex package (v0.3.0)}