Glmulti and Liner Mixed Models

glmulti and liner mixed models

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

glmulti , lmer fit (linear mixed models) and gls fit models(lme package)

This works for me with sessionInfo() as follows:

R Under development (unstable) (2014-09-17 r66626)
Platform: i686-pc-linux-gnu (32-bit)

other attached packages:
[1] glmulti_1.0.7 rJava_0.9-6   lme4_1.1-8    Rcpp_0.11.2   Matrix_1.1-4 

loaded via a namespace (and not attached):
[1] compiler_3.2.0  grid_3.2.0      lattice_0.20-29 MASS_7.3-34    
[5] minqa_1.2.3     nlme_3.1-117    nloptr_1.0.4    splines_3.2.0  
[9] tools_3.2.0  

With the following code:

library("lme4")
library("glmulti")
dd <- read.table("SO_glmulti.dat",header=TRUE)
m1 <- lmer(Yeild~ (TShann+Alt+Slope+CPT+MAT+MARF)^2+
               (1|Blocks)+(1|Composition),
           data=dd)

Note I get a warning about predictor scaling here -- probably harmless

From previous SO question:

setMethod('getfit', 'merMod', function(object, ...) {
    summ <- coef(summary(object))
    summ1 <- summ[,1:2,drop=FALSE]
    ## 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 is the old version of glmulti -- quick and dirty but depends on deparsing the formula.

lmer.glmulti<-function(formula,data,random="",...) {
    lmer(paste(deparse(formula),random),data=data,
         REML=FALSE,...)
}

Harder to grok but more robust:

lmer.glmulti<-function(formula,data,random="",...) {
    newf <- formula
    newf[[3]] <- substitute(f+r,
                            list(f=newf[[3]],
                                 r=reformulate(random)[[2]]))
    lmer(newf,data=data,
         REML=FALSE,...)
}

This is what I ended up with:

glmulti_lmm <- glmulti(formula(m1,fixed.only=TRUE),
                     random="+(1|Blocks)+(1|Composition)",
                       data=dd,method="g",
                       deltaM=0.5, 
                       fitfunc=lmer.glmulti,
                       intercept=TRUE,marginality=FALSE,level=2)

I initially tried the default method="h", gave up after 2650 models. My first run with method="g" got to a fairly stable IC after 50 generations, but the mean IC kept going down slowly, so I got impatient and boosted deltaM to 0.5.

For the first run I got IC=1651.69761603866 with the model Yeild~1+CPT+CPT:TShann+CPT:Alt+MAT:Alt+MARF:CPT

On the second run (with deltaM increased), I got a little bit luckier (IC=1649.61044009369, Yeild~1+TShann+CPT+CPT:Alt+MAT:CPT+MARF:CPT). (I don't know whether there's a way to set seed/ensure reproducibility with glmulti). The model claimed to converge after 120 generations.

coef(glmulti_lmm) worked fine for me. The bottom of the output (highest-weighted variables) were:

                  Estimate Uncond. variance Nb models   Importance
[... skip ...]
CPT:MARF      1.334119e-03     5.491836e-07        11 0.7875438701
CPT:MAT       4.051261e-02     5.215084e-04        18 0.7995790960
TShann        1.260082e+00     1.650145e+00        35 0.8111493166
CPT          -9.923303e-01     2.764638e-01        74 0.9600979205
Alt:CPT      -2.465917e-04     7.937155e-09        72 0.9765910742
(Intercept)   3.754893e+01     5.988814e+01       100 1.0000000000

By the way, you might be interested in "Ecologists overestimate the importance of predictor variables in model averaging: a plea for cautious interpretations", Galipaud et al. http://dx.doi.org/10.1111/2041-210X.12251

Related Topics