42

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 http://glmm.wikidot.com/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 = 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 mean it isn't going to tell us much anyway?

Kerry
  • 1,129
  • 3
  • 14
  • 20
  • 1
    Regarding obtaining an exact 0 variance of a random effect see https://stats.stackexchange.com/questions/115090. – amoeba Jan 21 '18 at 21:00

2 Answers2

27

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; ie. your model is degenerate.

As you correctly identify yourself: most probably, yes; ID as a random effect is unnecessary. 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 it is on the boundaries of the parameters 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 (eg. 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 of the 3 but actually gives you the best idea really about what is 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.

EdM
  • 57,766
  • 7
  • 66
  • 187
usεr11852
  • 33,608
  • 2
  • 75
  • 117
  • 7
    The RLRsim package is a really convenient way to test random effects using simulation-based likelihood ratio tests. – atrichornis Apr 15 '13 at 14:46
  • 1
    @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. – usεr11852 Apr 15 '13 at 15:08
  • 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! – atrichornis Apr 15 '13 at 18:30
  • 5
    The bootstrap on the boundary of the parameter space [has its own set of issues and problems leading to inconsistency](http://www.citeulike.org/user/ctacmo/article/1227040). The bootstrap is not a panacea, and should not be thrown in to the bag lightly assuming it will resolve everything. – StasK Apr 17 '13 at 13:34
  • @StasK : thank you for the reference; I will definitely have a look at it (and probably update my answer accordingly). I didn't mean to imply that parametric bootstrap is magical statistical elixir. Nevertheless I believe that realistically it is better than any other readily available method (AIC, anova) because at least it does not assume a specific underlying distribution for the statistic of interest. – usεr11852 Apr 17 '13 at 17:48
  • 3
    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. – StasK Apr 18 '13 at 15:16
  • I wasn't able to go through the paper properly yet. As you say it is somewhat subtle; especially given that the paper examples refer to means and not variances (as we do in for this problem) I don't find it easy to directly interpreter those findings in the current case. (I mean, the variances do define an F-distribution distribution themselves being all positive, which is not the case in the paper settings) I am looking in Jiming Jiang, 2007 Chapt.2 for something a bit more generic. (I don't say there is no issue, I just need to understand the issue better!) – usεr11852 Apr 18 '13 at 16:14
  • @atrichornis: It seems to me that SAS glosses over the same uncertainties. I'd be interested in hearing if it were otherwise. – russellpierce Jul 26 '13 at 01:19
  • @StasK: The `RLRsim` package can only test a single variance component. If we are interested in testing the random slope (not assuming the correlation between the random intercept and random slope is 0), is there any potential solution or references? – Randel May 25 '14 at 00:05
  • how does `anova()` help you with random effect? It only shows fixed effects. – Tomas Jul 09 '19 at 19:12
  • @Curious: It allows for model comparison as it works with the difference in deviance. e.g. See Bates (2010) [`lme4: Mixed-effects modeling with R`](http://webcom.upmf-grenoble.fr/LIP/Perso/DMuller/M2R/R_et_Mixed/documents/Bates-book.pdf) Sect. 2.2.4. – usεr11852 Jul 10 '19 at 08:17
3

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 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 that I'm currently using).

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.

VLC
  • 379
  • 2
  • 6
  • 1
    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`) – usεr11852 Apr 16 '13 at 23:41
  • @user11852 Please, see my update and let us know if this approach works also with `lme4_0.99999911-1`. – VLC Apr 17 '13 at 08:27
  • 1
    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. – usεr11852 Apr 17 '13 at 12:17
  • @user11852 Doesn't it work or won't it work? In other words, did you try it with `lme4_0.99999911-1`? In Galwey's book (pp. 202-203) the same approach (with identical results) is used also with another software (GenStat). – VLC Apr 17 '13 at 14:59
  • 4
    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. – usεr11852 Apr 17 '13 at 17:08
  • @user11852 I'm happy with your answer, I'm still wondering why Galwey suggests this approach. Maybe one day he'll read this post and share his opinion. – VLC Apr 17 '13 at 18:10
  • Indeed. Assuming that `lmer()` reflects Bates, Maechler and Bolker's insights then probably they thought the same originally but later changed their mind about it. :D (That's why earlier versions allow this but newer ones don't) Panta rhei. – usεr11852 Apr 17 '13 at 18:17
  • @usεr11852 See an answer that was just posted. It should probably better be a comment here, I will flag accordingly. – amoeba Jan 21 '18 at 21:06
  • @amoeba: I agree. – usεr11852 Jan 21 '18 at 23:55
  • 1
    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: https://github.com/lme4/lme4/issues/411. – Robin Beaumont Jan 21 '18 at 20:10
  • Even though you can operationally bypass the check for a constant effect, it is conceptually unnecessary since you can simply go back to a plain vanilla linear model and there are tools to do that available. Can you just NOT include the random effect in the lmer call? Failing that, can anova compare a simple lm output to a lmer output? – W7GVR Oct 30 '20 at 16:20