Confidence bands for model averaged predictions of GLMMs I use R with the MuMIn package for Multimodel inference.
My global Model is
M.global <- glmer(y ~ c1 + I(c1^2) + c2 + c3 + c4 + (1|subject),data, family =poisson(link="log))

After that I use
M.dredge <- dredge(M.global, m.max=3, subset = dc(cpl,I(cpl^2)))

To calculate all possible models from the global model and rank them by AICc. m.max limits the number of predictor to 3 and dc is a dependency chain, so the function takes the squared c1 only if unsquared c1 is included. m.max is used since the total sample size is 30, and according to Harrels 10:1 rule of thumb there is amximum of 3 predictors to avoid overfitting.
The I use 
M.avg <- model.avg(M.dredge, subset= delta < 4, fit=TRUE))
to get the model average across a subset of candidate models with a maximum difference of 4 AICc points to the "best" model having the lowest AICc. fit = TRUE refits every model in the subset.
Now I can use M.avg to get the averaged (full and conditional) coefficients for all predictors included in the subset as well as for the confdidence intervals:
Coef <- coef(M.avg, full=TRUE)
CI <- confint(M.avg, full=TRUE)

Now my question: Using poisson GLMMs and model averaging, how can I calculate model predictions and according confidence intervals for these predictions, to be used in a graph showing the regression line and confidence bands for that line?
Here they fill in the coefficients in the Model formula by hand and calculate point estimates for a given value (0, and 0.25) for the predictor f. Afterwards they repeat the procedure with the lower and the upper value of the confidence interval for the coefficients. However, I repeated that and I got something like this:

 A: Hi thanks to lots of hints by Kamil Bartón, i ended up with the following procedure:
M.global <- glmer(y ~ c1 + I(c1^2) + c2 + c3 + c4 + (1|subject),data, family =poisson(link="log))
M.dredge <- dredge(M.global, m.max=3, subset = dc(cpl,I(cpl^2)))
M.avg <- model.avg(M.dredge, subset= delta < 4, fit=TRUE))

newdata <- expand.grid(c1= with(data, seq(min(c1), max(c1), length.out=50)),
                       c2 = mean(data$c2),
                       c3 = mean(data$c3),
                       c4 = mean(data$c4))
# No random effects variable is needed!

# Calculate Predicitons with stnadard errors, and confidence intervals subsequently

pred.se <- predict(M.avg, type="response", se.fit=TRUE, re.form=NA, full=T, newdata)
  newdata$fit <- pred.se$fit
  newdata$SE <- pred.se$se
  newdata$upr=newdata$fit+1.96*newdata$SE
  newdata$lwr=newdata$fit-1.96*newdata$SE

# full=T defines the use of  full model-averaged coefficients see [here][1] 
# (full=F would not work either, due to mixed models); type= "response" 
# calculates predictions on the response scale (i.e. uses the inverse link
# function, exp() in poisson); re.form = NA avoids the use of random factors
# (inclusion is not possible with se.fit=TRUE).

Note: when using glmmadmb(family="nbinom2") the use of type="response" is not possible, but type="link" is. You could take the exp() of fit, SE, lwr, and upr, afterwards in the newdata frame.
The plots can be produced by:
p1 <- ggplot(newdata, aes(x=c1, y=fit)) +
    geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .25) +
    geom_line(size = 2) +
    labs(x = "c1 continous variable (restandardized)", y = "y response variable on original scale")
  print(p1)

Don't foget to re-standadize z transformed predictor variables in the final newdata frame. i.e.:
x <- attr(data, 'scaled:scale')
scale <- data.frame(t(x))
x <- attr(data, 'scaled:center')
center <- data.frame(t(x))

newdata$c1 <- newdata$c1 * scale$c1 + center$c1
newdata$c2 <- newdata$c2 * scale$c2 + center$c2
newdata$c3 <- newdata$c3 * scale$c3 + center$c3
newdata$c4 <- newdata$c4 * scale$c4 + center$c4

Now I am just wondering, why calculation of confidence intervals from single general linear mixed models (GLMMs) are so intensively discussed. Often you hear that you need complicated bootstrap methods or MarkovChainMonteCarlo (MCMC) algorithms to produce confidence intervals. Very, very complex for non-mathematicians. And even the function confint() takes much longer to calculate profile confidence intervals from a single mixed models than does from type=averaging model objects. Probably this is a Question of Accuracy and the se.fit = TRUE method is the least accurate? Furthermore i don't know if a correction method for Confidence intervals is included in the procedure presented here.
