Attempt to produce Huber-White standard errors for Poisson glmer() fits

robust.R

# Score residuals
modelRes <- function(model) residuals(model, "pearson")^2 / residuals(model, "working")
# Working residuals
# modelRes <- function(model) residuals(model, "working")

estfunGlmm <- function(x, ...)
{
  xmat <- model.matrix(x)
  xmat <- naresid(na.action(x), xmat)
  if(any(alias <- is.na(fixef(x)))) xmat <- xmat[, !alias, drop = FALSE]
  wres <- as.vector(modelRes(x)) * weights(x, "working")
  rval <- wres * xmat
  return(rval)
}

breadGlmm <- function(x, ...)
{
  sx <- summary(x)
  wres <- as.vector(modelRes(x)) * weights(x, "working")
  return(as.matrix(vcov(x) * as.vector(length(wres))))
}

sandwichGlmm <- function(x, Meat, ...)
{
  Bread <- breadGlmm(x)
  n <- nrow(estfunGlmm(x))
  return(1/n * (Bread %*% Meat %*% Bread))
}

robustSEglmm <- function(model, cluster){
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- length(fixef(model))
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj <- apply(estfunGlmm(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwichGlmm(model, Meat = crossprod(uj)/N)
  return(as.matrix(sqrt(diag(rcse.cov))))
}