CODEWITHSUNDEEP
rlme4glmulti
R.Timmers
R.Timmers

Reputation: 3

Strange glmulti results: Why are interaction variables from the candidate model dropped/not included?

I have been using glmulti to obtain model averaged estimates and relative importance values for my variables of interest. In running glmulti I specified a candidate model for which all variables and interactions were included based on a priori knowledge (see code below).

After running the glmutli model I studied the results by using the functions summary() and weightable(). There seem to be a number of strange things going on with the results which I do not understand.

First of all, when I run my candidate model with lme4 glmer() function I obtain an AIC value of 2086. In the glmulti output this candidate model (with exactly the same formula) has a lower AIC value (2107), as a result of which it appears at position 8 out of 26 in the list of all potential models (as obtained through the weigtable() function).

What seems to be causing this problem is that the logArea:Habitat interaction is dropped from the candidate model, despite level=2 being specified. The function summary(output_new@objects[[8]]) provides a different formula (without the logArea:Habitat interaction variable) compared to the formula provided through weightable(). This explains why the candidate model AIC value is not the same as obtained through lme4, but I do not understand why the interaction variables logArea:Habitat is missing from the formula. The same is happening for other possible models. It seems that for all models with 2 or more interactions, one interaction is dropped.

Does anyone have an explanation for what is going on? Any help would be much appreciated!

Best, Robert

Note: I have created a subset of my data (https://drive.google.com/open?id=1rc0Gkp7TPdnhW6Bw87FskL5SSNp21qxl) and simplified the candidate model by removing variables in order to decrease model run time. (The problem remains the same)

     newdat <- Data_ommited2[, c("Presabs","logBodymass", "logIsolation", "Matrix", "logArea", "Protection","Migration", "Habitat", "Guild", "Study","Species", "SpeciesStudy")] 

     glmer.glmulti <- function (formula, data, random, ...) {
     glmer(paste(deparse(formula), random), data = data, family=binomial(link="logit"),contrasts=list(Matrix=contr.sum, Habitat=contr.treatment, Protection=contr.treatment, Guild=contr.sum),glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000)))
     }

   output_new <- glmulti(y = Presabs ~  Matrix  + logArea*Protection + logArea*Habitat,
    data = sampledata,
    random = '+(1|Study)+(1|Species)+(1|SpeciesStudy)',
    family = binomial,
    method = 'h',
    level=2,
    marginality=TRUE,
    crit = 'aic',
    fitfunc = glmer.glmulti,
    confsetsize = 26)

    print(output_new)

    summary(output_new)

    weightable(output_new)

Upvotes: 0

Views: 422

Answers (1)

R.Timmers
R.Timmers

Reputation: 3

I found a post (https://stats.stackexchange.com/questions/341356/glmulti-package-in-r-reporting-incorrect-aicc-values) of someone who encountered the same problem and it appears that the problem was caused by this line of code: 

glmer.glmulti <- function (formula, data, random, ...) {
  glmer(paste(deparse(formula), random), data = data, family=binomial(link="logit"))
}

By changing this part of the code into the following the problem was solved:

glmer.glmulti<-function(formula,data,random,...) {
  newf <- formula
  newf[[3]] <- substitute(f+r,
                          list(f=newf[[3]],
                               r=reformulate(random)[[2]]))
  glmer(newf,data=data,
        family=binomial(link="logit"))
}

Upvotes: 0

Related Questions