title: How to add confidence intervals for a zero-inflated negative binomial function in ggplot2 R?
body: <p>I have a data set that is modelled as zero-inflated negative binomial with mixed effects. I want to get the confidence intervals from the model predictions and plot the models mean and confidence intervals. I have attempted to plot the model mean. Can someone let me know if this is the correct way to do it? I do not know how to plot the confidence intervals from the models onto ggplot2. I would like to plot the predicted mean of the data along with it's confidence interval. My basic attempt at the graph has been made below in the code.</p>
  
  <pre><code>library(pscl)
  library(lmtest)
  
  df &lt;- data.frame(
        fertilizer = c("N","N","N","N","N","N","N","N","N","N","N","N","P","P","P","P","P","P","P","P","P","P","P","P","N","N","N","N","N","N","N","N","N","N","N","N","P","P","P","P","P","P","P","P","P","P","P","P"), 
        level = c("low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","high","low","low","high","low"), 
        repro = c(0,90,2,4,0,80,1,90,2,33,56,0,99,100,66,80,1,0,2,33,0,0,1,2,90,5,2,2,5,8,0,1,90,2,4,66,0,0,0,0,1,2,90,5,2,5,8,55)
      )    
  
      model &lt;- formula(repro ∼ fertilizer + level | fertilizer * level)
      modelzinb &lt;- zeroinfl(model, dist = "negbin", link = "logit",data =df)
      summary(modelzinb)
  
      df$predict &lt;- predict(modelzinb)
  
      ggplot(df, aes(x=fertilizer, y=predict, color = fertilizer)) + theme_bw() +  stat_summary(aes(color = fertilizer),fun.y = mean, geom = "point", size = 4, position = position_dodge(0.1)) +
        scale_x_discrete("Fertlizer") +
        facet_wrap(.~level)  
  </code></pre>