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?


2 Answers 2


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.

  • 10
    $\begingroup$ The RLRsim package is a really convenient way to test random effects using simulation-based likelihood ratio tests. $\endgroup$ Apr 15, 2013 at 14:46
  • 2
    $\begingroup$ @atrichornis: +1. Interesting package; I did not know about it. I just had a look at its code, quite straightforward I might say. I wish they incorporate it (or something like that) to lme4 especially now that mcmcsamp() is broken and people are left with only their own ad-hoc bootstrap implementations to get some decent p-values etc. out. $\endgroup$
    – usεr11852
    Apr 15, 2013 at 15:08
  • $\begingroup$ True, mixed models are not straightforward in R. Lots of approximations and workarounds... Although I gather SAS etc just gloss over some of the same uncertainties? Ben Bolker is a coauthor on both packages, he may have his reasons for not incorporating it. Probably time! $\endgroup$ Apr 15, 2013 at 18:30
  • 6
    $\begingroup$ The bootstrap on the boundary of the parameter space has its own set of issues and problems leading to inconsistency. The bootstrap is not a panacea, and should not be thrown in to the bag lightly assuming it will resolve everything. $\endgroup$
    – StasK
    Apr 17, 2013 at 13:34
  • 3
    $\begingroup$ Do take a look, the argument is very subtle. As far as I can recall, it boils down to the fact that you are doing the bootstrap from a distribution that is different from the null; and given the non-standard distributions obtained on the boundary, the regularity conditions are violated, and the bootstrap distribution does not converge to the target. I think non-parametric bootstrap can still be constructed here by taking out the group means of residuals. However, with violation of independence of observations between groups, another layer of complications may arise. $\endgroup$
    – StasK
    Apr 18, 2013 at 15:16

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)


@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).

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)))

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

anova(fm1, fm0)


  Data: dataset
  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.

  • 1
    $\begingroup$ Have you tested your idea? Please prove me wrong but I think your idea won't even execute. If the IDconst is the same for all your data then you don't have any grouping. You need a grouping factor to have at least one sampled level and the way you set the model up it has none. I could maybe believe the rationale of using a "random grouping" but that's a different ball-game all together. Test your approach with some dummy data. I strongly believe that with your proposed setup lmer() won't run. (I use lme4_0.99999911-1) $\endgroup$
    – usεr11852
    Apr 16, 2013 at 23:41
  • $\begingroup$ @user11852 Please, see my update and let us know if this approach works also with lme4_0.99999911-1. $\endgroup$
    – VLC
    Apr 17, 2013 at 8:27
  • 1
    $\begingroup$ Nope, it won't work. And I said, "it shouldn't" because conceptually you are not having a mixed model to start with. (Sorry I might sounded to aggressive in my previous comment.) I guess what you want to try to achieve is set up the $Z$ matrix to be a unitary vector but unless you find a way to explicitly do that in R (ie. write your own deviance function) you are out of luck. $\endgroup$
    – usεr11852
    Apr 17, 2013 at 12:17
  • 4
    $\begingroup$ Yes, I did what you suggest; it will not work/compute. Error in lFormula(formula = value ~ fac1 + (1 | idconst), data = dataset) : grouping factors must have at least 1 sampled level. And as I said, conceptually it is wrong. It is not a matter of tricking the software to give some numbers out, it is a matter if what you say it reasonable. You don't have a second mixed model to compare with if in that model the random effect is by construction a constant. You might as well exclude it and try a linear model instead. $\endgroup$
    – usεr11852
    Apr 17, 2013 at 17:08
  • 1
    $\begingroup$ Update concerting defining a single group random variable in lme4. This can be done if you set the option: control=lmerControl(check.nlev.gtr.1="ignore"). Ben Bolker mentions it here: github.com/lme4/lme4/issues/411. $\endgroup$ Jan 21, 2018 at 20:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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