tools for using glmulti with lme4/glmmTMB model fits

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")