-
-
Save bbolker/4ae3496c0ddf99ea2009a22b94aecbe5 to your computer and use it in GitHub Desktop.
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") | |
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)".
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)