How to plot mixed-effects model estimates in ggplot2 in R?

I have a 2x2x2 factorial design with one random effect. The data (dat) is as follows:

  colour  size  level   marbles set
  Blue    Large Low     80      1
  Blue    Large High    9       2
  Blue    Small Low     91      1
  Blue    Small High    2       1 
  White   Large Low     80      2
  White   Large High    9       1
  White   Small Low     91      2
  White   Small High    2       1

I want to plot two models:

mod1 <- lmer(marbles ~ colour + size + level + colour:size + colour:level + size:level + (1|set), data = dat)

mod2 <- lmer(marbles ~ colour + size + level +(1|set), data = dat)

I usually use the following code to do my plots:

pd <- position_dodge(0.82)
  ggplot(dat, aes(x=colour, y=marbles, fill = level)) + theme_bw() + 
  stat_summary(geom="bar", fun.y=mean, position = "dodge") +  
  stat_summary(geom="errorbar", fun.data=mean_cl_boot, position = pd)+
  + facet_grid(~size)

I'm unsure on how to replace the terms with coefficients from the model estimates. Any ideas on how can I plot the estimates of the final model in gpplot2? It would be helpful if anyone can suggest a easy way to print the model estimates too

In addition, is there anyway that I can get ggplot2 to display bars on top of the graphs showing interactions that are significant?