1

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")
Kara
  • 6,115
  • 16
  • 50
  • 57
matmar
  • 318
  • 2
  • 13
  • 1
    The first error I get is: `pippo<-as.data.frame(cbind(vy2,va,vb,random_effect)) Error in cbind(vy2, va, vb, random_effect) : object 'vy2' not found` – IRTFM Dec 06 '13 at 17:54

2 Answers2

1

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

Vincent Calcagno
  • 26
  • 1
0

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.

IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • Thanks for the reply. I'm sorry for the error in the code. My purpose is to obtain the model parameter coefficients averaged between all the models in the glmulti.fit. This is what the coef.glmulti does. Unfortunately what you showed doesn't help me since it only shows the parameters of the best model. I am pretty convinced that the problem is in the structure of the lmer.fit object that differs from lm.fit. So the function cannot extract the coefficients of each models using the getfit() function. I will follow your suggestion writing to the authors of the package. – matmar Dec 06 '13 at 20:36
  • How does one "average coefficients" across mutiple models with different covariates? – IRTFM Dec 06 '13 at 20:58
  • There are two versions of model averaging: 1) The coefficients beta of the variable x is averaged over all models in which the variable x appers; 2) The alterantive is to consider that variable x is in every model, it is just that in some models the corresponding beta value is set to zero, rather than being considered unknown. There is a broad literature about model selection and multi-model inference, if you are interested I suggest you to have a look at Burnham and Anderson 2002. – matmar Dec 07 '13 at 17:50
  • OK. I suppose one could do such a maneuver; and I will do some reading; it's just that my experience in the biomedical domain is that coefficients are generally attenuated in "larger" models relative to their univariate counterparts, positive correlations of risk factors being the norm. I admit that I have avoided the use of automatic procedures in general, having been taught to prefer domain knowledge as a guide when available. I'm just an old fogey I suppose. – IRTFM Dec 07 '13 at 18:16