Multiple Regression Especification in one Plots

I'd like to make a plot like in R

But, I have no idea how to do it. Can some of you help me? The plot shows the coefficients and Confidence intervals for several regressions and each specification is identify by several characteristics like the use of controls, the type of polynomial approximation and bandwidth method (in a RD design). The paper where I saw the graph can be found here https://smontanoc.github.io/assets/files/BMM_SignalingSkills.pdf and it is the Appendix Figure D.2.

This looks like a customized type of upset plot. Without knowing your starting point, it's very difficult to give concrete advice. However, we may show an example of how you could go about creating such a plot.

Suppose we had two different models that estimated a continuous outcome based on a categorical variable:

library(tidyverse)
library(patchwork)

mod1 <- lm(price ~ cut, data = within(diamonds, cut <- as.character(cut)))
mod2 <- glm(price ~ cut, data = within(diamonds, cut <- as.character(cut)),
            family = Gamma)

We just need to create a data frame of predictions from both models:

cuts <- c("Fair", "Good", "Very Good", "Premium", "Ideal")

preds1 <- predict(mod1, se.fit = TRUE, newdata = data.frame(cut = cuts))
preds2 <- predict(mod2, se.fit = TRUE, newdata = data.frame(cut = cuts))

df <- tibble(model = rep(c("linear", "Gamma"), each = 5),
           cut = factor(rep(cuts, 2), cuts),
           fit = c(preds1$fit, 1/preds2$fit),
           upper = c(preds1$fit + qnorm(0.975) * preds1$se.fit,
                     1/(preds2$fit + qnorm(0.975) * preds2$se.fit)),
           lower = c(preds1$fit + qnorm(0.025) * preds1$se.fit,
                     1/(preds2$fit + qnorm(0.025) * preds2$se.fit)),
           combo = fct_inorder(interaction(model, cut)))

Now we create three plots: one to show our estimates, one to indicate our categorical variables, and one to indicate the model used:

p1 <- ggplot(df, aes(combo, fit)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1) +
  theme_minimal() +
  labs(y = "Price", x = NULL) +
  theme(axis.text.x = element_blank())

p2 <- ggplot(df, aes(x = combo, y = cut)) +
  geom_point(color = "gray", size = 4,
             data = expand.grid(combo = df$combo, cut = df$cut)) +
  geom_point(size = 4) +
  theme_void() +
  ggtitle("Cut") +
  theme(axis.text.y = element_text(hjust = 1),
        plot.title.position = "plot")

p3 <- ggplot(df, aes(x = combo, y = model)) +
  geom_point(color = "gray", size = 4,
             data = data.frame(combo = rep(df$combo, 2),
                               model = rep(df$model, each = 2))) +
  geom_point(size = 4) +
  theme_void() +
  ggtitle("Model") +
  theme(axis.text.y = element_text(hjust = 1),
        plot.title.position = "plot")

Finally, we combine the plots using patchwork:

p1 / plot_spacer() /  p2 / plot_spacer() / p3 + 
  plot_layout(heights = c(20, 1, 5, 1, 2))