How can I test whether a random effect is significant?

Question

I am trying to understand when to use a random effect and when it is unnecessary. Ive been told a rule of thumb is if you have 4 or more groups/individuals which I do (15 individual moose). Some of those moose were experimented on 2 or 3 times for a total of 29 trials. I want to know if they behave differently when they are in higher risk landscapes than not. So, I thought I would set the individual as a random effect. However, I am now being told that there is no need to include the individual as a random effect because there is not a lot of variation in their response. What I can't figure out is how to test if there really is something being accounted for when setting individual as a random effect. Maybe an initial question is: What test/diagnostic can I do to figure out if Individual is a good explanatory variable and should it be a fixed effect - qq plots? histograms? scatter plots? And what would I look for in those patterns.

I ran the model with the individual as a random effect and without, but then I read Ben Bolker's GLMM FAQ where they state:

do not compare lmer models with the corresponding lm fits, or glmer/glm; the log-likelihoods are not commensurate (i.e., they include different additive terms)

And here I assume this means you can't compare between a model with random effect or without. But I wouldn't really know what I should compare between them anyway.

In my model with the random effect I also was trying to look at the output to see what kind of evidence or significance the RE has:

lmer(Velocity ~ D.CPC.min + FD.CPC + (1|ID), REML=FALSE, family=gaussian, data=tv)

Linear mixed model fit by maximum likelihood 
Formula: Velocity ~ D.CPC.min + FD.CPC + (1 | ID) 
   Data: tv 
    AIC    BIC logLik deviance REMLdev
 -13.92 -7.087  11.96   -23.92   15.39
Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 0.00000  0.00000 
 Residual             0.02566  0.16019 
Number of obs: 29, groups: ID, 15

Fixed effects:
              Estimate Std. Error t value
(Intercept)  3.287e-01  5.070e-02   6.483
D.CPC.min   -1.539e-03  3.546e-04  -4.341
FD.CPC       1.153e-04  1.789e-05   6.446

Correlation of Fixed Effects:
          (Intr) D.CPC.
D.CPC.min -0.010       
FD.CPC    -0.724 -0.437

You see that my variance and SD from the individual ID as the random effect equals 0. How is that possible? What does 0 mean? Is that right? Then my friend who said "since there is no variation, using ID as random effect is unnecessary" is correct? So, then would I use it as a fixed effect? But wouldn't the fact that there is so little variation suggest it isn't going to tell us much anyway?

Answer 1

The estimate ID's variance = 0, indicates that the level of between-group variability is not sufficient to warrant incorporating random effects in the model; i.e., your model is degenerate.

As you correctly identify yourself: most probably, yes; ID as a random effect is unnecessary. A few things spring to mind to test this assumption:

1. You could compare (using REML = F always) the AIC (or your favourite IC in general) between a model with and without random effects and see how this goes.
2. You would look at the anova() output of the two models.
3. You could do a parametric bootstrap using the posterior distribution defined by your original model.

Mind you, choices 1 & 2 have an issue: you are checking for something that is on the boundaries of the parameter space so actually they are not technically sound. Having said that, I don't think you'll get wrong insights from them and a lot of people use them (e.g., Douglas Bates, one of lme4's developers, uses them in his book but clearly states this caveat about parameter values being tested on the boundary of the set of possible values).

Choice 3 is the most tedious but actually gives you the best idea about what is really going on. Some people are tempted to use non-parametric bootstrap also but I think that given the fact you are making parametric assumptions to start with you might as well use them.

Answer 2

I'm not sure that the approach I'm going to suggest is reasonable, so those who know more about this topic correct me if I'm wrong.

My proposal is to create an additional column in your data that has a constant value of 1:

IDconst <- factor(rep(1, each=length(tv$Velocity)))

Then, you can create a model that uses this column as your random effect:

fm1 <- lmer(Velocity ~ D.CPC.min + FD.CPC + (1|IDconst), 
            REML=FALSE, family=gaussian, data=tv)

At this point, you could compare (AIC) your original model with the random effect ID (let's call it fm0) with the new model that doesn't take into account ID (fm1) since IDconst is the same for all your data.

anova(fm0, fm1)

Update

@user11852 was asking for an example, because in his/her opinion the above approach won't even execute. On the contrary, I can show that the approach works (at least with lme4_0.999999-0).

set.seed(101)
dataset <- expand.grid(id=factor(seq_len(10)), fac1=factor(c("A", "B"),
    levels=c("A", "B")), trial=seq_len(10))
dataset$value <- rnorm(nrow(dataset), sd=0.5) +
  with(dataset, rnorm(length(levels(id)), sd=0.5)[id] +
  ifelse(fac1 == "B", 1.0, 0)) + rnorm(1, .5)
dataset$idconst <- factor(rep(1, each=length(dataset$value)))

library(lme4)
fm0 <- lmer(value ~ fac1 + (1|id), data=dataset)
fm1 <- lmer(value ~ fac1 + (1|idconst), data=dataset)

anova(fm1, fm0)

Output:

  Data: dataset
  Models:
  fm1: value ~ fac1 + (1 | idconst)
  fm0: value ~ fac1 + (1 | id)

      Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
  fm1  4 370.72 383.92 -181.36                      
  fm0  4 309.79 322.98 -150.89 60.936      0  < 2.2e-16 ***

According to this last test, we should keep the random effect since the fm0 model has the lowest AIC as well as BIC.

Update 2

By the way, this same approach is proposed by N. W. Galwey in Introduction to Mixed Modelling: Beyond Regression and Analysis of Variance on pages 213-214.