Skip to content

Instantly share code, notes, and snippets.

@bbolker
Last active March 26, 2024 00:54
Show Gist options
  • Save bbolker/4ae3496c0ddf99ea2009a22b94aecbe5 to your computer and use it in GitHub Desktop.
Save bbolker/4ae3496c0ddf99ea2009a22b94aecbe5 to your computer and use it in GitHub Desktop.
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")
@bbolker
Copy link
Author

bbolker commented Mar 3, 2024

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

@bbolker
Copy link
Author

bbolker commented Mar 3, 2024

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

@gudryan
Copy link

gudryan commented Mar 3, 2024

@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.

@bbolker
Copy link
Author

bbolker commented Mar 3, 2024

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

@gudryan
Copy link

gudryan commented Mar 4, 2024

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.

@Study-yang666
Copy link

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.

@bbolker
Copy link
Author

bbolker commented Mar 23, 2024

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)

@Study-yang666
Copy link

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)

@bbolker
Copy link
Author

bbolker commented Mar 25, 2024

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

@Study-yang666
Copy link

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment