glmmTMB_GH774.R

library(MuMIn)
library(lme4)
library(glmmTMB)
library(emmeans)

set.seed(2) ## simulate some data...
dd <- expand.grid(f1 = factor(1:3), f2 = factor(1:3),
                  g = factor(1:10), rep = 1:2)
Z <- model.matrix(~g-1, dd)
b <- rnorm(10)
X <- model.matrix(~f1 + f2, dd)
beta <- 1:5
eta <- X%*%beta + Z%*%b
dd$y <- rnbinom(nrow(dd), mu = exp(eta), size = 1)

options(na.action = "na.fail")
mod1 <- glmer.nb(y ~ f1 + f2 + (1|g), data = dd)
mod1_d <- dredge(mod1, rank = "AIC")
mod1_a <- model.avg(mod1_d, fit = TRUE, data = dd)
emmeans(mod1_a, list(pairwise ~ f1 | f2), data=dd, adjust="Tukey")

mod2 <- glmmTMB(y ~ f1 + f2 + (1|g), data = dd, family = nbinom2)
mod2_d <- dredge(mod2, rank = "AIC")
mod2_a <- model.avg(mod2_d, fit = TRUE, data = dd)
emmeans(mod2_a, list(pairwise ~ f1 | f2), data=dd, adjust="Tukey")

##
vcov(mod2_a, full = TRUE)
MuMIn:::vcov.averaging(mod2_a, full = TRUE)


## hacked version of vcov.averaging
va <- function (object, full = FALSE, ...) {
    object <- MuMIn:::upgrade_averaging_object(object)
    full <- MuMIn:::.checkFull(object, full)
    full <- as.logical(full)[1L]
    models <- attr(object, "modelList")
    if (is.null(models))
      stop("cannot calculate covariance matrix from ", "'averaging' object without component models")
    fix_mat <- function(x) {
      v <- as.matrix(x$cond)
      nn <- sprintf("cond(%s)", rownames(v))
      nn <- gsub("(Intercept)", "Int", nn)
      dimnames(v) <- list(nn, nn)
      v
    }
    vcovs <- lapply(lapply(models, vcov), fix_mat)
    names.all <- dimnames(object$coefArray)[[3L]]
    nvars <- length(names.all)
    nvarseq <- seq(nvars)
    wts <- Weights(object)
    wts <- wts/sum(wts)
    vcov0 <- matrix(if (full)
        0
    else NA_real_, nrow = nvars, ncol = nvars, dimnames = list(names.all,
        names.all))
    vcovs2 <- lapply(vcovs, function(v) {
        i <- match(MuMIn:::fixCoefNames(dimnames(v)[[1L]]), names.all)
        vcov0[i, i] <- v
        return(vcov0)
    })
    b1 <- object$coefArray[, 1L, ]
    if (full)
        b1[is.na(b1)] <- 0
    avgb <- object$coefficients[2L - full, ]
    res <- sapply(nvarseq, function(c1) sapply(nvarseq, function(c2) {
        weighted.mean(sapply(vcovs2, "[", c1, c2) + (b1[, c1] -
            avgb[c1]) * (b1[, c2] - avgb[c2]), wts, na.rm = TRUE)
    }))
    dimnames(res) <- list(names.all, names.all)
    return(res)
}

va(mod2_a, full = TRUE)

emmeans(mod2_a, list(pairwise ~ f1 | f2), data=dd, adjust="Tukey",
        vcov. = va(mod2_a, full = TRUE))
## still fails