title: Plotting multiple model specifications on same ggplot
body: <p>I'm working with a linear model and I'd like to show how robust the effect is to different specifications using a single point range plot of differences in means across conditions. This is as far as I've gotten (MWE below).</p>
  <p>I have three important dummy treatment indicators, plus five covariates.</p>
  <p>What I want to do now is overlay estimates of the three treatments in this plot with plots for five additional models, each with a different covariate included, then add a legend/shapes/colors to distinguish them. I assume I could <code>group_by()</code> and <code>do()</code> five separate models, but confidence interval names get replaced, and I'm not sure how to get ggplot to read more than one model (especially in tidyverse, which is foreign to me).</p>
  <p>I wasn't able to figure it out or find any existing threads that deal quite with a problem like this. Can this be done? Thanks for any tips!</p>
  <p>MWE with example data:</p>
  <pre><code>treatment1 = rep(seq(0, 1, 1), 300)
  treatment2 = sample(seq(from = 0, to = 1, by = 1), size = 300, replace = TRUE)
  treatment3 = rep(seq(0, 1, 1), each=300)
  response = rnorm(n = 300, mean = 3, sd = 1)
  cov1 = rnorm(n = 300, mean = 0, sd = 1)*response
  cov2 = rnorm(n = 300, mean = 0, sd = 1)/response
  cov3 = rnorm(n = 300, mean = 0, sd = 1)-response
  cov4 = rnorm(n = 300, mean = 0, sd = 1)+response
  cov5 = rnorm(n = 300, mean = 0, sd = 1)*log(response)
  df &lt;- data.frame(treatment1,treatment2,treatment3,
             response,cov1,cov2,cov3,cov4,cov5)
  
  model &lt;- df %&gt;% group_by(treatment1, treatment2,treatment3) %&gt;% 
    do(data.frame(tidy(lm(response ~ treatment1*treatment2*treatment3, data = .),conf.int=T, conf.level = 0.95 )))
  
  facet.labs &lt;- c(&quot;T1=0&quot;, &quot;T1=1&quot;)
  names(facet.labs) &lt;- c(&quot;0&quot;, &quot;1&quot;)
  
  
  model$treatment3 &lt;-factor(model$treatment3, labels = c(&quot;T3=1&quot;,&quot;T3=0&quot;))
  model$treatment3 &lt;-factor(model$treatment3, levels = c(&quot;T3=1&quot;,&quot;T3=0&quot;))
  
  ggplot(model, aes(x=estimate, y=treatment2, shape = treatment3)) +
    geom_pointrange(position = position_dodge(width = 1), aes(xmin=conf.low, xmax=conf.high), size=.75) +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(), 
          panel.background = element_blank(),
          axis.line = element_line(colour = &quot;black&quot;),
          panel.border = element_rect(colour = &quot;black&quot;, fill=NA, size=.8)) +
    scale_y_continuous(name =&quot;&quot;, breaks = c(.1,1.22), labels=c(&quot;T2=0&quot;,&quot;T2=1&quot;)) +
  geom_vline(xintercept=0, linetype=&quot;dotted&quot;) +
    labs(title=&quot;&quot;) + xlab(&quot;&quot;)  +
    labs(shape=&quot;&quot;)+ theme(axis.ticks = element_blank()) +
    theme(axis.text.y = element_text(angle = 90, vjust = 0.5, hjust=1)) +#,col=&quot;Treatment 1&quot;)+
    guides(fill = guide_legend(override.aes = list(linetype = 0,fill=NA)))+
    facet_wrap(~treatment1, labeller = labeller(treatment1 = facet.labs)) +
    guides(shape = guide_legend(reverse=T))
  
  </code></pre>