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