tools for using glmulti with lme4/glmmTMB model fits

glmulti_mixed.R

library(lme4)
library(glmmTMB)
library(glmulti)
set.seed(101)

## updated 3 March 2024; take random effects as formulas so that e.g. `ar1()` terms work

## FIXME: better data so we don't get convergence weirdness etc.?
## A random vector of count data
vy1 <- round(runif(100, min=1,max=20)*round(runif(100,min=1,max=20)))

## Predictors
va <- runif(100, min=1, max=100)
vb <- runif(100, min=1, max=100)
random_effect <- factor(rep(c(1,2,3,4),each=25))
tt <- factor(rep(1:25, 4))
pippo <- data.frame(vy1,va,vb,random_effect, tt)
form_glmulti <- vy1~va*vb

## The wrapper function for linear mixed-models
mm_glmulti <- function(formula, data, random=NULL, FUN = lme4::lmer,
                       extra_args = NULL, ...) {
    ## FIXME: check/throw error if REs not parenthesized?
    af <- function(x,y) if (is.null(y)) x else glmmTMB::addForm(x,y)
    ff <- af(formula, random)
    arglist <- c(list(ff, data = data), list(...), extra_args)
    do.call(FUN, arglist)
}

lmer.glmulti <- function(...) mm_glmulti(..., extra_args = list(REML = FALSE))
glmer.glmulti <- function(...) mm_glmulti(..., FUN = lme4::glmer)
glmmTMB.glmulti <- function(...) mm_glmulti(..., FUN = glmmTMB::glmmTMB,
                                             extra_args = list(REML = FALSE))


## I'm not sure why I set `df` to 10000 but I assume it's just to get Gaussian
##  confidence intervals (and that `Inf` doesn't work for some reason)
setMethod('getfit', 'merMod',
          function(object, ...) {
              summ <- coef(summary(object))
              summ1 <- summ[, c("Estimate", "Std. Error"), drop = FALSE]
              ## FIXME: update for lmerTest?
              cbind(summ1, df=rep(10000, length(fixef(object))))
          }
          )

setMethod('getfit', 'glmmTMB',
          function(object, ...) {
              ## only handles conditional
              ## warn if disp, zi detected??
              summ <- coef(summary(object))$cond
              summ1 <- summ[, c("Estimate", "Std. Error"), drop = FALSE]
              cbind(summ1, df=rep(10000, length(fixef(object)$cond)))
          }
          )

## sigh, can't collapse this to a single function that passes
## fitfunc because of glmulti() eval()/call() stuff
##  (Error in match.fun(fitfunction) : object 'FUN' not found)
glmulti_lmm <- glmulti(form_glmulti,
                        random= ~(1|random_effect),
                        data=pippo, method="h", 
                        fitfunc=lmer.glmulti,
                        intercept=TRUE,marginality=FALSE,level=2)

glmulti_glmm <- glmulti(form_glmulti,
                        random= ~(1|random_effect),
                        family = poisson,
                        data=pippo, method="h", 
                        fitfunc=glmer.glmulti,
                        intercept=TRUE,marginality=FALSE,level=2)

glmulti_glmmTMB <- glmulti(form_glmulti,
                        random= ~(1|random_effect),
                        data=pippo, method="h", 
                        fitfunc=glmmTMB.glmulti,
                        intercept=TRUE, marginality=FALSE,level=2)

glmulti_glmmTMB_nbinom <- glmulti(form_glmulti,
                                  random= ~(1|random_effect),
                                  family = "nbinom2",
                                  data=pippo, method="h", 
                                  fitfunc=glmmTMB.glmulti,
                                  intercept=TRUE, marginality=FALSE,level=2)

glmulti_glmmTMB_ar1A <- glmulti(form_glmulti,
                           random=~(ar1(tt-1|random_effect)),
                           data=pippo, method="h", 
                           fitfunc=glmmTMB.glmulti,
                           intercept=TRUE, marginality=TRUE,level=2)

coef.glmulti(glmulti_lmm,
             select="all",
             varweighting="Johnson",
             icmethod="Burnham")

coef.glmulti(glmulti_glmmTMB,
             select="all",
             varweighting="Johnson",
             icmethod="Burnham")


coef.glmulti(glmulti_glmmTMB_nbinom,
             select="all",
             varweighting="Johnson",
             icmethod="Burnham")

Yes, I got the same error when using "random = ~(1 | pair.fate$Year)" or "random= ~(1|Year)+(1|PairID)".