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

Hi everyone,
Do these wrappers work just for linear models?
I'm trying to use your wrapper glmmTMB in negative binomial models and it is not working.. receive the error:

Error in terms.formula(formula, data = data) :
Invalid model formula in ExtractVars

What did you try? Can we have a reproducible example? With the code in the original gist, this seems to work:

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)
coef.glmulti(glmulti_glmmTMB_nbinom,
             select="all",
             varweighting="Johnson",
             icmethod="Burnham")

@jukop, I finally got around to updating this so that your ar1() example works.

@bbolker , thank you very much for your reply. Unfortunately, I have no reproducible example.. I tried a lot but could not to suggest any data set that could reveal the same problem.

It seems to be related exactly to the wrapper function and not to how it is included in the glmulti function.
I trying to apply as follow:

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

If I include all other possible distribution, such as poison:

glmer.glmulti<-function(formula,data,random="",...){glmer(paste(deparse(formula),random),data=data,...)}

it works perfectly. Just negative binomial from glmmTMB package did work. What could cause these error?

Error in terms.formula(formula, data = data) :
Invalid model formula in ExtractVars

Now, trying your suggestion I get the same error message again.

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

Thanks or posting the gist for using glmulti with glmmTMB. Two notes I encountered:

1. I changed af <- function(x,y) if (is.null(y)) x else glmmTMB::addForm(x,y) to af <- function(x,y) if (is.null(y)) x else reformulas::addForm(x,y) to avoid this error:

Error: 'addForm' is not an exported object from 'namespace:glmmTMB'

2. When I run

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

I get the following error --

in method for ‘getfit’ with signature ‘"glmmTMB"’: no definition for class “glmmTMB”

So I tried this and it ran:

setOldClass("glmmTMB")
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)))
          }
)