1
$\begingroup$

First of all, I am a beginner in statistics, and I have been reading some questions and comments here regarding this matter but I am still a bit lost.

My problem: I am constructing GLMM's in order to assess habitat selection. For now I am doing some exploratory analysis to each one of my variables. I am checking AIC (and AICc, which has been the same for AIC) and R squared (both using MuMIn R package), and overdispersion, using RVAideMemoire R package. I am a beginner in statistics, and I have been reading some questions and comments here regarding this matter but I am still a bit lost.

Regaring the dataset:

X is the dependent variable, it is binomial (1 and 0); 1 represents "used" locations, obtained through telemetry techniques, and 0 represents "available", which are random points collected across each study area. Variable A represents study areas and Variable B represents interactions between individuals (represented in codes, for instance: M01-M02; M01-M02-F01, and so on), and each group of individuals is unique for each study area. Variable 1 represents a given type of land cover.

str(data1)
'data.frame':   32670 obs. of  8 variables:
 $ VarA: Factor w/ 5 levels
 $ VarB: Factor w/ 51 levels 
 $ X      : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Var1     : num  -0.201 -0.201 -0.201 -0.201 -0.201 ...
 $ Var2     : num  1.383 -0.973 -0.748 1.611 -0.973 ...
 $ Var3     : num  -0.0985 -0.0985 -0.0985 -0.0985 -0.0985 ...
 $ Var4     : num  -0.482 -0.482 0.942 -0.482 -0.482 ...
 $ Var5     : num  -0.502 0.293 0.813 -0.783 2.41 ...

For each variable, I am constructing a GLMM using two different random effects, because I don't know which one is best/most appropriate for my data. For instance, I have these models for the same variable:

    lm1 <- glmer(x~Var1+(1|VarA), data=data1, family = binomial(link="logit"))
    lm1_1 <- glmer(x~Var1+(1|VarA/VarB), data=data1, family = binomial(link="logit"))
    lm1_2 <- glmer(x~Var1+(1|VarA)+(1|VarA/VarB), data=data1, family = binomial(link="logit"))

For each model, I check AIC, R squared and overdispersion, and I don't know how to interpretate my results. For instance:

lm1: AIC-19908.9; overdispersion-0.609;r2m and r2c almost zero (4.06e-07)

lm1_1: AIC-8806.5; overdispersion-0.261;r2m-0.0016 and r2c-0.7923

For lm1_2 I get:

Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.00174688 (tol = 0.001, component 1)

ANOVA for lm1 and lm1_1:

> anova(lm1, lm1_1)
Data: data1
Models:
lm1: Used ~ Var1 + (1 | VarA)
lm1_1: Used ~ Var1 + (1 | VarA/VarB)
        Df     AIC     BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
lm1    3 19908.9 19934.1 -9951.4  19902.9                            
lm1_1  4  8806.5  8840.1 -4399.3   8798.5 11104      1  < 2.2e-16 ***

This pattern repeats itself across other variables.

1) How do I interpret this? Because, if I focus only on AIC, lm1 is the most appropriate, and also if I consider overdispersion. But R2 is terrible. If I focus on R2, it is only good for lm1_1 and only for the nested random effect, which would indicate that the nested random effect is "correct". Am I right?

2) Are there any other good approaches for GLMM?

Improve this question
$\endgroup$
  • $\begingroup$ Can you provide more details about your data - for instance what are VarA and VarB ? Normally the physical/biological situation would dictate whether or not to include higher levels of grouping. Trying to select a model based on a combination of AIC and psuedo $R^2$ is often futile as they are both difficult to apply to GLMMs and they are assessing different things. Also, note that (1|VarA/VarB) is equivalent to ` (1 | VarA)+(1 | VarA:VarB) ` which will fit a random intercept for both VarA and VarB with VarB within VarA so the (1|VarA) in your formula for lm1_1 is redundant. $\endgroup$ – Robert Long Jul 3 '16 at 13:44
  • $\begingroup$ I am sorry, I had the formula for lm1_1 wrong, now it is update. I only considered for lm1_1 "(1|VarA/VarB)", that is why I was confused. When I used both, I got the message "model failed to converge" (updated as well). $\endgroup$ – mtao Jul 3 '16 at 14:44
  • $\begingroup$ Please can you also explain what X is and also how the data were collected. Is each individual unique within each area ? Also please provide information about the dataset (total observations, total individuals, total areas) perhaps include a str() of the data. $\endgroup$ – Robert Long Jul 3 '16 at 14:52
  • $\begingroup$ It is updated now!! $\endgroup$ – mtao Jul 3 '16 at 15:08
  • $\begingroup$ Why are variables Var2-Var5 not in the model ? Also, please post the output from glmer() for lm1 and lm1_1 along with the output from anova(lm1,lm1_1) $\endgroup$ – Robert Long Jul 3 '16 at 15:27
2
$\begingroup$

You have only 5 levels of VarA, which is too low to use it as a random effect. This can lead to incorrect estimates and numerical problems. Since its estimated variance from the output that you provided in chat was 1.226e-08 this is another reason not to include it.

You have measured 4 other variables which are of importance in your study, apart from Var1 so these should be included as covariates, subject to considerations of mediation and collinearity. You might consider including one or more of these as random slopes.

As for model selection, I think we have established that lm1_1 can be discarded, but if/when you have competing models it is not unusual to find contradiction from pseudo-$R^2$ and AIC. Since you have a large enough dataset you can do cross-validation, where you build your models on one subset of the data and assess their predictive abilities on a different subset.

Finally, overdispersion appears to be minimal and not a cause for concern.

Improve this answer
$\endgroup$
  • $\begingroup$ I conducted some analyses with the same data and I had what seem to be weird results. I constructed a model with all variables (fixed/random effects), and then used "dredge()" to select the "best" models, according to AICc delta and weight. When I use VarB (ID), I obtain very low values for weight, but VarB as a random effect has variance 13.02 (SD 3.6),which seems significant. But if I conduct the same test but with a "glm()" instead of "glmer", thus without the random effect, weight values are higher. What does it mean? Should I not use a random effect at all? Or is not related? $\endgroup$ – mtao Jul 9 '16 at 16:55
  • $\begingroup$ @Teresa, please ask a new question about this, giving as much information as possible. $\endgroup$ – Robert Long Jul 9 '16 at 17:37
  • $\begingroup$ done! here: stats.stackexchange.com/questions/222949/… $\endgroup$ – mtao Jul 9 '16 at 18:30
2
$\begingroup$

Without any information on your data, it's difficult to help you properly. One thing that struck me is that in you 2nd model (lm1_1) you have 2 times the same random term: 1 time VarA is not nested and the 2nd time VarB is nested in VarA. I believe you'll have to decide which one to keep, is VarB nested in VarA (then only keep "+(1|VarA/VarB)" or is it not (then only keep "+(1|VarA)"). I can't think of a reason why you would have those 2 in your model.

If you are a beginner, you probably should spend some time reading and understanding how to build GLMM. I suggest you to follow the steps described in the book from Zuur et al. (2009) Mixed Effects Models and Extensions in Ecology with R to understand how do build a GLMM. I read it 2 times and it really helped me to understand GLMMs.

Improve this answer
$\endgroup$
  • 1
    $\begingroup$ I don't think it makes any difference, the (1|VarA) in lm1_1 is simply redundant. +1 for the Zuur reference ! $\endgroup$ – Robert Long Jul 3 '16 at 13:47
  • $\begingroup$ @RobertLong Good to know that (1|VarA) does not make a difference and is simply redundant. However, I would say that it's always good to keep your codes as simple as possible for clarity. $\endgroup$ – Mud Warrior Jul 3 '16 at 14:16
  • $\begingroup$ I updated my post, I wrote wrongly lm1_1, sorry, it only includes "(1|VarA/VarB)". When I include both, I get the message "model failed to converge". $\endgroup$ – mtao Jul 3 '16 at 14:46
  • $\begingroup$ Well, that's probably because of this redundant (1|VarA) term. But really, you should look into Zuur's book, GLMMs are complex beasts and I don't know how I would have navigated through this without the explanations provided in this book. I would also add that never in Zuur et al. book I remember them saying to base your best model on R2. If you don't have to add more explanatory variables, I would say that lm1_1 is the best so far (lowest AIC), your overdispersion is fine (lower than 1) and if individuals are nested within sites than it seems fine to me. $\endgroup$ – Mud Warrior Jul 3 '16 at 15:02
  • $\begingroup$ Is it ok if the overdispersion value is low, for instance, 0.2? Is there an interval of "best" values for this measure? $\endgroup$ – mtao Jul 3 '16 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.