title: Confidence Intervals for negative binomial GLMM: uncertainty in the variance parameters
body: <p>I'd like to calculate the 95 confidence Intervals for prediction in negative binomial GLMMs, using <code>lme4</code> package. But, a cause of no option for computing standard errors of predictions by default I try to:</p>
  <pre><code>#Package
  library(lme4)
  
  # Create a data set
  set.seed(101)
  dd &lt;- expand.grid(f1 = factor(1:3),
                    f2 = LETTERS[1:2], g=1:9, rep=1:15,
                    f3 = rnorm(10,25),
            KEEP.OUT.ATTRS=FALSE)
  mu &lt;- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2)))          
  dd$y &lt;- rnbinom(nrow(dd), mu = mu, size = 0.5)
  str(dd)
  
  # Create a negative binomial generalized mixed model
  m.nb &lt;- glmer.nb(y ~ f1 + f2 + poly(f3,2) + (1|g), data=dd, verbose=TRUE)
  m.nb
  
  #Compute the confidence interval for model predictions
  mm&lt;-model.matrix(~f1+f2+poly(f3,2),dd)
  dd$y_est&lt;-mm%*%fixef(m.nb) #predicted values
  pvar &lt;- diag(mm %*% tcrossprod(vcov(m.nb),mm))
  CIupper &lt;- dd$y_est+sqrt(pvar) * 1.96 
  CIlower &lt;- dd$y_est-sqrt(pvar) * 1.96
  </code></pre>
  <p>But I am not so sure if this is an efficient method for incorporating uncertainty in the variance parameters, especially using a <code>poly()</code> variable. Please, any help with this?</p>