Interpreting random effect variance in glmer

Question

I'm revising a paper on pollination, where the data are binomially distributed (fruit matures or does not). So I used glmer with one random effect (individual plant) and one fixed effect (treatment). A reviewer wants to know whether plant had an effect on fruit set -- but I'm having trouble interpreting the glmer results.

I've read around the web and it seems there can be issues with directly comparing glm and glmer models, so I'm not doing that. I figured the most straightforward way to answer the question would be to compare the random effect variance (1.449, below) to the total variance, or the variance explained by treatment. But how do I calculate these other variances? They don't seem to be included in the output below. I read something about residual variances not being included for binomial glmer -- how do I interpret the relative importance of the random effect?

> summary(exclusionM_stem)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(Fruit_1, Fruit_0) ~ Treatment + (1 | PlantID)

     AIC      BIC   logLik deviance df.resid 
   125.9    131.5    -59.0    117.9       26 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0793 -0.8021 -0.0603  0.6544  1.9216 

Random effects:
 Groups  Name        Variance Std.Dev.
 PlantID (Intercept) 1.449    1.204   
Number of obs: 30, groups:  PlantID, 10

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.5480     0.4623  -1.185   0.2359   
TreatmentD   -1.1838     0.3811  -3.106   0.0019 **
TreatmentN   -0.3555     0.3313  -1.073   0.2832   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr) TrtmnD
TreatmentD -0.338       
TreatmentN -0.399  0.509

score 13 · Accepted Answer · answered May 23 '15 at 02:24

13

While getting an analogue of "proportion variance explained by each effect" is in principle possible for GLMMs, there are several complicating factors (which levels of the model do you consider "total variance", and how do you quantify the sampling variation due to the lowest-level [Binomial in this case] sampling distribution)? Nakagawa and Schielzeth (doi:10.1111/j.2041-210x.2012.00261.x) present a general approach to calculating R^2 (proportion of total variance explained) for (G)LMMs which has gotten pretty popular in ecology; Xu et al 2003 take a similar approach. In principle this approach could probably be extended to consider proportion of variance explained by different terms [but note that the 'proportion of variance' of all the terms in the model considered in this way would probably not add up to 100% -- it could be either more or less].

However, if your reviewer isn't hung up on statistical details and would be satisfied with a more heuristic explanation of "importance", you could point out that the estimated among-plant standard deviation is 1.20, very close to the magnitude of the largest treatment effect (-1.18); this means that the plants vary quite a bit, relative to the magnitude of the treatment effects (e.g., the 95% range of the plant effects is approximately $4\sigma$, from $-1.96 \sigma$ to $+1.96\sigma$).

Visually:

enter image description here

answered May 23 '15 at 02:24

Ben Bolker

34,308
2
93
126

+1, I'm intrigued by your (I infer favorable) mention of $R^2$ for non-linear models. What is your take on the discussion here: [Which pseudo-R2 measure is the one to report for logistic regression (Cox & Snell or Nagelkerke)?](http://stats.stackexchange.com/q/3559/7290) – gung - Reinstate Monica May 23 '15 at 02:35
1

I was just saying that I don't think it's a crazy or necessarily ill-posed question. But both hierarchical structure and GLM-like models open cans of worms that make it harder to pick an answer. I usually don't bother, but I can see why people would want to try find numbers that measured goodness of fit, or relative importance of terms in a model, in a reasonable way. – Ben Bolker May 23 '15 at 02:47
That's reasonable. Btw, what do you think of my suggestion that w/ 10 plants, 3 treatments & N=30, the OP could fit a model using both as fixed effects? I don't necessarily think this would be the right final model, of course, but it does strike me as a potentially permissible way to test if there is variation amongst plants, & to put both variables on similar footing for comparison. – gung - Reinstate Monica May 23 '15 at 02:59
seems reasonable to me. – Ben Bolker May 23 '15 at 03:30
I fit a model with both Treatment and Plant as fixed effects as gung suggested, and the Plant term had very high p-value (p = 0.3). Does this seem odd given that, as you say, "the estimated among-plant standard deviation is 1.20, very close to the magnitude of the largest treatment effect (-1.18)"? Why would it show up as insignificant in an ANOVA with 2 fixed effects? – jwb4 May 23 '15 at 03:47

score 3 · Answer 2 · answered May 22 '15 at 19:30

What you want is to test if the variance of PlantID is $0$. However, this is a weird test to try to run, because null value is at the boundary of the permissible space. Such tests are still run, but a lot of people are very uncomfortable with them.

In your case, you have multiple measures per plant, so one quick and dirty approach is to run a model with PlantID as a fixed effect, and test that effect.

score 1 · Answer 3 · edited Apr 13 '17 at 12:44

The simple answer to your reviewer is, "Yes." If he is asking you to test whether the variance of the random effect is significantly different from 0, you have a couple options. Note though that many smart people are uncomfortable with testing if variances of random effects are different from 0.

Simplest is a likelihood ratio test, though not recommended by most. They are very conservative when testing at the boundaries (i.e. your are testing against a variance of 0 which is as low as it can go). There is a rule of thumb out there that the p value is about twice what it really is.

The method recommended most places is a parametric bootstrap. You can use bootMer from the lme4 package. Make sure that your set the REML parameter of your lmer function to FALSE, otherwise your variance will be greater than 0 100% of the time (or close to it... actually it'll probably be greater than 0 nearly 100% of the time anyway).

Some tips and further resources:

http://glmm.wikidot.com/faq (find the How can I test whether a random effect is significant? heading)

lmer() parametric bootstrap testing for fixed effects

http://www.r-bloggers.com/using-bootmer-to-do-model-comparison-in-r/

Thanks for this lucid (and prompt!) guide to model comparison. But how would I interpret the "magnitude" of the effect of the random variable? i.e., how would I compare the variance explained by my random variable to the variance explained by the fixed variable (treatment)? I guess I don't see how this is gleaned from the results of the bootstrapped LRT test. — jwb4, May 23 '15 at 00:28

score 1 · Answer 4 · edited Apr 13 '17 at 12:44

1

In Multiple-Sample Cochran's Q Test , they use anova to compare the results of the two models (one without random effects and one with random effects).

Jairo Rocha University of the Balearic Islands

edited Apr 13 '17 at 12:44

Community

1

answered Mar 10 '17 at 12:34

Jairo Rocha

11
1

Interpreting random effect variance in glmer

4 Answers4

Linked