CODEWITHSUNDEEP
r
matmar
matmar

Reputation: 348

glmulti and liner mixed models

I have a problem with the glmulti package and linear mixed models. When I try to estimate the model-averaged coefficients with the coff.glmulti function I get this error:

Error in data.frame(..., check.names = FALSE) :arguments imply differing number of rows: 1, 0

I did a little bit of debugging and I found that the problem starts in the highlighted line of the coeff.glmulti function:

         if (length(object@adi) >= 1) 
             for (j in 1:length(object@adi)) {
               cak[[length(names(cak)) + 1]] = object@adi[[j]]

               names(cak)[length(names(cak))] = names(object@adi)[j]
             }
         modf = eval(cak)

         coffee = c(coffee, list(modf))
     }
 }

 if (length(coffee) == 1) {
     warning("Only one candidate: standard conditional inference was performed.")

     **return(coef(coffee[[1]]))**
 }

After, when it tries to apply getfit on the coffee object it fails. I think that the error is due to a different structure of the lmer.fir object respect to lm or other type of models object.

I'm pasting a minimum repeatable example to facilitate who wants to help me:

#Add the required package
library(lme4)
library(glmulti)

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

# The wrapper function for linear mixed-models
lmer.glmulti<-function(formula,data,random="",...){
  lmer(paste(deparse(formula),random),data=data,REML=F,...)
}
# The wrapper function for linear models
lm.glmulti<-function(formula,data,...){
  lm(paste(deparse(formula)),data=data,...)
}

# Multi selection for lmer
glmulti_lmm<-glmulti(form_glmulti,random="+(1|random_effect)",data=pippo,method="h", 
                     fitfunc=lmer.glmulti, intercept=TRUE,marginality=FALSE,level=2)
# Model selection for lm
glmulti_lm<-glmulti(form_glmulti,data=pippo,method="h",fitfunc=lm.glmulti,intercept=TRUE, 
                        marginality=FALSE, level=2)

# Coeffs estimation lmer #Here the error
coef.glmulti(glmulti_lmm,select="all",varweighting="Johnson",icmethod="Burnham")
#Coeffs estimation lm #With lm everything is ok
coef(glmulti_lm,varweighting="Johnson",icmethod="Burnham",select="all")

Upvotes: 2

Views: 3648

Answers (2)

Vincent Calcagno
Vincent Calcagno

Reputation: 36

If you are using one of the latest versions of lme4, the getfit() function I recommended is no longer adapted. Indeed, the lme4 package maintainers have made quite a lot of changes in their package: the class of objects is now "merMod", where it was "mer", and a few other things.

Then the getfit function must be slightly adjusted in order to interface glmulti with the new lme4 structure. Here a a getfit definition that works with the latest builds of lme4 for Ubuntu 12.04, as of yesterday:

setMethod('getfit', 'merMod', function(object, ...) {
summ=summary(object)$coef
summ1=summ[,1:2]
if (length(dimnames(summ)[[1]])==1) {
    summ1=matrix(summ1, nr=1, dimnames=list(c((Intercept)"),c("Estimate","Std.    Error")))
}
cbind(summ1, df=rep(10000,length(fixef(object))))
})

This should fix the issue. [see also my website http://vcalcagnoresearch.wordpress.com/package-glmulti/] Regards

Upvotes: 2

IRTFM
IRTFM

Reputation: 263451

After changing the token name from vy2 to vy1 the first call to glmulti no longer throws an error but this one does:

coef.glmulti(glmulti_lmm,select="all",varweighting="Johnson",icmethod="Burnham")

I do not think that your highlighted line is the source of the error since if it were the warning would have been issued and furthermore the object has 8 models. I think it is just below that point where this line appears:

coke = lapply(coffee, getfit) # since that is step three in the traceback()

Looking at the content of glmulti_lmm we see that the objects slot has 8 models:

> summary(glmulti_lmm)$bestmodel
[1] "vy1 ~ 1"

> glmulti_lmm@objects[[1]]
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: vy1 ~ 1 + (1 | random_effect) 
   Data: data 
      AIC       BIC    logLik  deviance 
1183.9965 1191.8120 -588.9983 1177.9965 
Random effects:
 Groups        Name        Std.Dev.
 random_effect (Intercept)  0.00   
 Residual                  87.45   
Number of obs: 100, groups: random_effect, 4
Fixed Effects:
(Intercept)  
      105.5  

> coef( glmulti_lmm@objects[[1]])
$random_effect
  (Intercept)
1      105.48
2      105.48
3      105.48
4      105.48

attr(,"class")
[1] "coef.mer"

You didn't say what your goal was but perhaps this shows how to inspect such objects. To me this is suspicious for a bug and you may want to send a memo to the package maintainer.

Upvotes: 1

Related Questions