3
$\begingroup$

I am developing GLMM's in order to assess habitat selection (using GLMMs' coeficients to construct Resource selection functions). I have (telemetry) data from 5 study areas, and each area has a different number of individuals monitored.

To develop GLMM's, the dependend variable is binary (1-used locations; 0-available locations), and I have a initial set of 14 continuous variables (8 land cover variables; 2 distance variables, to artificial areas and water sources; 4 topographic variables): a buffer was placed around each location and the area of each land cover within that buffer was accounted for; distances were measured from each point to the nearest feature, and topographic variables were obtained using DEM rasters. I tested for correlation using Spearman's Rank, so not all 14 were used in the GLMMs. All variables were transformed using z-score.

As random effect, I used individual ID (In another question ("GLMM: relationship between AIC, R squared and overdispersion?"), it became clear that using study areas as random effect was not useful nor correct).

I constructed a GLMM with 9 variables (not correlated) and a random effect, then used "dredge()" function and "model.avg(dredge)" to sort models by AIC values. This was the result (only models of AICc lower than 2 represented):

[1]Call:
model.avg(object = dredge.m1.1)

Component model call: 
glmer(formula = Used ~ <512 unique rhs>, data = All_SA_Used_RP_Area_z, family = 
     binomial(link = "logit"))

Component models: 
          df   logLik    AICc  delta weight
123578     8 -4309.94 8635.89   0.00   0.14
1235789    9 -4309.22 8636.44   0.55   0.10
123789     8 -4310.52 8637.04   1.14   0.08
1235678    9 -4309.75 8637.50   1.61   0.06
12378      7 -4311.78 8637.57   1.67   0.06
1234578    9 -4309.79 8637.58   1.69   0.06

Variables 1 and 2 represent the distance variables; from 3 to 8 land cover variables, and 9 is a topographic variable. Weights seem to be very low, even if I average all those models as it seems to be common when delta values are low. Even with this weights, I constructed GLMMs for each of the combinations, and the results were simmilar for all 6 combinations. Here are the results for the first one (GLMM + overdispersion + r-squared):

Random effects:
 Groups    Name        Variance Std.Dev.
 ID.CODE_1 (Intercept) 13.02    3.608   
Number of obs: 32670, groups:  ID.CODE_1, 55

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.54891    0.51174  -1.073 0.283433    
3       -0.22232    0.04059  -5.478 4.31e-08 ***
5       -0.05433    0.02837  -1.915 0.055460 .  
7       -0.13108    0.02825  -4.640 3.49e-06 ***
8       -0.15864    0.08670  -1.830 0.067287 .  
1         0.28438    0.02853   9.968  < 2e-16 ***
2         0.11531    0.03021   3.817 0.000135 ***     
Residual deviance: 0.256           
r.squaredGLMM():
       R2m        R2c 
0.01063077 0.80039950 

This is what I get from this analysis:

1) Variance and SD of the random effect seems fine (definitely better than the "0" I got when using Study Areas as random effect);

2) Estimate values make sense from what I know of the species and the knowledge I have of the study areas;

3) Overdispersion values seem good, and R-squared values don't seem very good (at least when considering only fixed effects) but, as I read in several places, AIC and r-squared are not always in agreement.

4) Weight values seem very low. Does it mean the models are not good?

Then what I did was construct a GLM ("glm()"), so no random effect was used. I used the same set of variables used in [1], and here are the results (only models of AICc lower than 2 represented):

[2] Call:
model.avg(object = dredge.glm_m1.1)

Component model call: 
glm(formula = Used ~ <512 unique rhs>, family = binomial(link = "logit"), data = 
     All_SA_Used_RP_Area_z)

Component models: 
          df   logLik     AICc   delta weight
12345678   9 -9251.85 18521.70    0.00   0.52
123456789 10 -9251.77 18523.54    1.84   0.21
1345678    8 -9253.84 18523.69    1.99   0.19

In this case, weight values are higher.

Does this mean that it is better not to use a random effect? (I am not sure I can compare GLMM with GLM results, correct me if I am doing wrong assumptions)

$\endgroup$
  • $\begingroup$ Can you update the question to include more detail about how the study is designed, including what the 14 independent variables are, and how they relate to each other. I know you said that it is telemetry data, but a little more information would be useful. $\endgroup$ – Robert Long Jul 9 '16 at 20:44
  • $\begingroup$ Done!! Do you need any more information? $\endgroup$ – mtao Jul 9 '16 at 22:53
  • $\begingroup$ Are you more interested in prediction (predicting habitat selection for new telemetry data), or inference (understanding habitat selection using telemetry data) ? $\endgroup$ – Robert Long Jul 10 '16 at 10:02
  • $\begingroup$ Kind of both. Fist I want to understand how these individuals selected the current habitat, and try to understand why. Then, with the coeficcients from GLMM I can build a RSF and project it in a given area, and there I can see which areas have high or low probability of finding the species. $\endgroup$ – mtao Jul 10 '16 at 10:31
2
$\begingroup$

I would strongly advise you to avoid automated model selection procedures such as dredge() (even the function name makes me shiver). There may be some merit in these when you are primarily concerned about prediction for future data, but even in this case it is strongly recommended to use some form of cross-validation, where you can build your model on a training dataset and then assess it's predictive capability on another dataset. If you build a model based on AIC with your whole dataset, while it may predict your current dataset well, there is a good chance it will perform poorly on new data. When your goal is mainly inference, the best way forward is to use theory and common sense to build your model. Unless you have a huge number of variables, then I think that theory and common sense is also a better method for prediction too.

A good starting point is to draw a path diagram to hypothesize the associations, and directions of causality, according to theory. This can allow you to build a model avoiding common problems such as over-adjustment for confounding, and including variables that should not be present in the model (for example, if they lie on the causal path between an exposure and the outcome). Although the current theory may not be not be well-developed (developing the theory is presumably one of your goals) a path diagram may help you rule out some possible models. DAGitty is a very user-friendly web-based graphical tool that can assist with this.

You have measured 14 variables, presumably choosing these for some good reasons. Excluding some solely on the basis of high bivariate correlation without understanding the relations between them is dangerous. If they are essentially measuring the same thing, or one is derived implicitly (or explicitly) from another, this may be valid, but if one is a cause of the other, then you need to think carefully about which one to exclude. A path diagram will be very useful for this.

You have repeated measures on individual subjects, therefore a priori it is a good idea to account for this by using random intercepts for subject because measurements on the same individual are likely to be more similar to those on another individual. I would strongly urge you not to rely on p-values in general, but even more so in this case. If you really want to test the significance of the random intercept then a bootstrap estimate is perhaps the most robust way to proceed.

As for your question about low weights, I don't know what these weights represent exactly, but I assume that they must all add up to 1 so if there are a lot of models to choose from then obviously there will be many models with small weights, and if the "best" ones are very similar to each other then by necessity, the "best" ones will have small weights. Note that the difference in AICc between the "best" and the worst (well, the 6th best, since there are only 6 shown) is just 1.69, which is telling you that there is very little difference between any of these models.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you very much for your help. So, what you are saying is that I should build my on models, based on the information I have about the species and predictions, and then still use "model.avg()" to access AICc, delta and weight parameters? But: when I do the "dredge()", I am using a full model with all variables (not correlated), and then dredge() gives me all possible combinations, so the combinations I construct will be for sure on that list, right? $\endgroup$ – mtao Jul 10 '16 at 11:10
  • $\begingroup$ No, I'm saying that you should try to avoid using these automated selection methods altogether unless you have a big enough dataset that you can split it and cross-validate. Even then, it is better to rely more on substantive knowledge than automated procedures. By all means, use AIC to choose between 2-3 models that you have arrived at via theory, but don't let AIC drive the whole model selection process. $\endgroup$ – Robert Long Jul 10 '16 at 13:29
  • $\begingroup$ Ok. What about the "weight" values? Do you have any idea of why are they so low? What can cause this? $\endgroup$ – mtao Jul 10 '16 at 13:45
  • $\begingroup$ @Teresa I just updated the answer $\endgroup$ – Robert Long Jul 10 '16 at 14:14
  • $\begingroup$ Thank you very much! Indeed I don't have enough data to perform cross-validation. So, even for predictions, I will construct several models with different sets of variables, based on my knowledge, and then run "model.avg()" again, let's see how the weights behave this time! $\endgroup$ – mtao Jul 10 '16 at 14:36
1
$\begingroup$

Akaike weights only provide information about the set of models from which they are calculated, so in your example, you can't really learn anything from comparing weights for the glmm set of models to weights for the glm set of models (i.e. Akaike weights won't tell you whether the random effect is appropriate).

There's a good discussion on testing random effects here. It does specifically state "do not compare lmer models with the corresponding lm fits, or glmer/glm; the log-likelihoods are not commensurate". That said, there is a worked glmer example by Ben Bolker here that does explicitly compare log-likelihoods and AICc values between glm and glmer models. If you do go this route, it will be the AICc values that you want to compare between glm and glmer, not the delta or weight values (which are only meaningful within a set).

The fact that the variance of your random effect is high (and sd relatively low) suggests to me that you should retain the random effect. Likewise, comparing the conditional and marginal R2 values (0.8 vs. 0.01) suggests that the random effect is explaining a lot of the variation in your response.

Update

It doesn't make sense to ask whether the glmer weight values are "too low". An Akaike weight is the probability that a model is the 'best', given the data and the set of models under consideration. If the top models in a given set all have low and similar Akaike weights, it just means that no one model in that set stands out as being much better than the other models in that set... it doesn't mean that the top models in your set are bad models in an absolute sense. It doesn't tell you either way. Akaike weights are only useful for relative comparison.

To assess goodness of fit in an absolute sense, you could look at R2 values. Your R2c of 80% suggests to me that that particular model is very good, though that's mostly driven by the random effect (i.e. individuals varied a lot in habitat use).

Since your aim is to get fixed effect coefficients for use in a RSF, and since no single model stands out as the 'best', model averaging would be a good approach. If mod is your glmer object, you can get model-averaged coefficients to use in your RSF as follows:

mod_dredge <- dredge(mod)
mod_avg <- model.avg(mod_dredge)
coefTable(mod_avg)

Since your fixed effects have relatively low explanatory power (R2m = 1%), an RSF based on those coefficients may not have much predictive power. As for what you can do to improve this model, it's tough to say... it could just be a biological reality that your study animals vary a lot in habitat use (thus the important random effect), but overall habitat use by these animals is simply not strongly related to the set of specific fixed effects you have examined.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for your answer! But don't you think the weight values (for the glmer) are too low? How can I justify this, or what can I do to improve the model? $\endgroup$ – mtao Jul 10 '16 at 8:30
  • $\begingroup$ @Teresa, see update $\endgroup$ – pbee Jul 11 '16 at 17:42
  • $\begingroup$ Thank you for your answer! I saw in an exercise about glmm where they average (model.avg()) the data twice. First they do that for the results using dredge(), which provides low weight values, and then they select the top models and average again, obtaining higher weight values. Which only shows that it is indeed a relative measure!! $\endgroup$ – mtao Jul 11 '16 at 22:18

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.