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

Hello,

I am trying to create a model with glmmTBM with ar1 covariance. Any ideas what need to be hacked into this code so that it would work? The function is printing me this:

Initialization...
TASK: Exhaustive screening of candidate set.
Error in ar1(tt - 1 | Org_ID) : could not find function "ar1"

Reproducible example please ... ?

Reproducible example please ... ?

Ok, I will do that and respond later with code. Thank you for your attention. :)

I simply added tt <- rep(1:25,4) into pippo data.frame and then modified later code to use ar1-structure. It appears that either code prints that ar1 function is not found or that 'ar1' is not an exported object from 'namespace:glmmTMB'.

I am not yet interested if the model would converge, just about technical stuff to make it work so that ar1 structure could be used.

I wonder if quasiquotation (https://adv-r.hadley.nz/quasiquotation.html) could be utilised to solve this problem. I have used it couple of times but I am not expert enough to make it work here. It also seems to have changed quite a lot since I used it last time three years ago.

# Author: bbolker, modified by jukop for discussion of use of ar1 structure
library(lme4)
library(glmmTMB)
library(glmulti)
set.seed(101)

## 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 <- rep(c(1,2,3,4),each=25)
tt <- rep(1:25,4)
pippo <- data.frame(vy1,va,vb,random_effect,tt=as.factor(tt))
form_glmulti <- vy1~va*vb

## The wrapper function for linear mixed-models
mm_glmulti <- function(formula, data, random="", FUN = lme4::lmer,
                       extra_args = NULL, ...) {
  ## FIXME: check/throw error if REs not parenthesized?
  ff <- glmmTMB::addForm(formula, substitute(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))

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)

# How to get ar1 to work with this code?
# ar1 structure (does not work)
glmulti_glmmTMB <- glmulti(form_glmulti,
                           random= ar1(tt-1|random_effect),
                           data=pippo, method="h", 
                           fitfunc=glmmTMB.glmulti,
                           intercept=TRUE, marginality=TRUE,level=2)
# ar1 structure second attempt
glmulti_glmmTMB <- glmulti(form_glmulti,
                           random= glmmTMB::ar1(tt-1|random_effect),
                           data=pippo, method="h", 
                           fitfunc=glmmTMB.glmulti,
                           intercept=TRUE, marginality=TRUE,level=2)
# so we need some way to tell computer to use ar1 structure from glmmTMB. Perhaps quasiquotation or something similar?

coef.glmulti(glmulti_glmmTMB,
             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)".