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

Can you please run traceback() after you get the error and post the results?

Can you show a little more of the code you are running? (e.g. what does your glmulti() call look like?)

Thank you @bbolker for your attention!

I solve it now and the error came from the way I was including the random effect. If I include:

random=(1|time)
it didn't works

if include
random=(1|dataset$time)
works perfectly

I was using the glmmTMB.glmulti wrapper function.

Hi! I followed the code above, but got error as: Error in glmulti(pair.rate.fit, random = ~(1 | pair.fate$Year), family = poisson, :
Improper call of glmulti. And setting "random=(1|dataset$time)" also did not work.

Can we have a reproducible example please? Are you able to run the example above exactly as specified? What version of glmulti are you using? (Version above tested with 1.0.8)

Yes, I can run the example, and it works. I use glmuti of version 1.0.8. Please see my code as following:

pair.rate.fit<-glmer(NumofFledglings ~
Mean.laying.date.median
+Difference.median
+Mean.laying.date.median:Difference.median
+Bodylengthf
+Tarsusf
+Billf
+Bodylengthm
+Tarsusm
+Billm
+Breeding.experience
+Age.f
+Age.m
+(1|Year)
+(1|PairID)
,data=pair.fledg
,family = poisson()
)

mm_glmulti <- function(formula, data, random=NULL, FUN = lme4::glmer,
extra_args = NULL, ...) {
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)
}

glmer.glmulti <- function(...) mm_glmulti(..., FUN = lme4::glmer)

setMethod('getfit', 'merMod',
function(object, ...) {
summ <- coef(summary(object))
summ1 <- summ[, c("Estimate", "Std. Error"), drop = FALSE]
cbind(summ1, df=rep(10000, length(fixef(object))))
}
)

glmulti_glmm <- glmulti(pair.rate.fit,
random= ~(1|Year)+(1|PairID) ,
family = poisson,
data=pair.fledg, method="h",
fitfunc=glmer.glmulti,
intercept=TRUE,marginality=FALSE,level=2)

I'll take a look if I get a chance. For what it's worth, the error that you report (which has "random = ~(1 | pair.fate$Year)")) doesn't seem to match the example you gave (which has "random= ~(1|Year)+(1|PairID)") ...

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