glmulti: assigning a predict function for glmer with two nested random variables

Question

I'm trying to use glmulti with glmer for model averaging and to get model averaged predictions. I've followed examples in the glmulti documentation ('Using glmulti with any type of statistical model, with examples', included with package) and updates provided on this website (glmulti and liner mixed models) and on the package maintainer's blog (https://vcalcagnoresearch.wordpress.com/package-glmulti/). I've managed to create a wrapper for the glmer function:

glmer.glmulti <- function (formula, data, family, random = "") {
glmer(paste(deparse(formula), random), data = data, family = binomial)
}

And I've managed to assign a getfit method for glmer so I can get the model averaged coefficients:

setMethod('getfit', 'merMod', function(object, ...) {
summ=summary(object)$coef
summ1=summ[,1:2]
if (length(dimnames(summ)[[1]])==1) {
summ1=matrix(summ1, nr=1, dimnames=list(c("(Intercept)"),c("Estimate","Std. Error")))
}
cbind(summ1, df=rep(10000,length(fixef(object))))
})

The next step is to assign a predict function for glmer. This is the example provided in the package documentation:

predict.mer=function(objectmer,random=random, newdata, withRandom=F,se.fit=F, ...){
if (missing(newdata) || is.null(newdata)) {
DesignMat <- model.matrix(objectmer) }
else {
DesignMat=model.matrix(delete.response(terms(objectmer)),newdata)
}
output=DesignMat %*% fixef(objectmer)
if(withRandom){
z=unlist(ranef(objectmer))
if (missing(newdata) || is.null(newdata)) {
Zt<- objectmer@Zt
} else {
Zt<-as(as.factor(newdata[,names(ranef(objectmer))]),"sparseMatrix")
}
output = as.matrix(output + t(Zt) %*% z)
}
if(se.fit){
pvar <- diag(DesignMat %*% tcrossprod(vcov(objectmer),DesignMat))
if(withRandom){
pvar <- pvar+ VarCorr(objectmer)[[1]]
}
output=list(fit=output,se.fit=sqrt(pvar))
}
return(output)
}

Then to get the the model averaged predictions (bab is a fitted glmulti object in the example):

> predict(bab, se.fit=T, withR=T)
Error in predict.merMod(coffee[[i]], se.fit = se.fit, ...) :
cannot calculate predictions with both standard errors and random effects

I've also tried:

> predict(bab, se.fit=T, withR=F)
Error in predict.merMod(coffee[[i]], se.fit = se.fit, ...) :
cannot calculate predictions with both standard errors and random effects

And:

> predict(bab, se.fit=F, withR=T)
Error in waou %*% t(matrix(unlist(preds), nrow = nbpo)) :
non-conformable arguments
In addition: There were 50 or more warnings (use warnings() to see the first 50)

I'm not quite sure what's wrong, although it may be something obvious. I gather that the lme4 package has been updated and changed since this example was written, so it may be something to do with that(?).

Another possibility is that the documentation says that this function will only handle one random variable. My model has two nested random variables: x ~ y + z + w + (1|u/v).

I need to a) get this working, b) update the function to handle two random variables. Any suggestions would be much appreciated.