36
$\begingroup$

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?

| cite | improve this question | |
$\endgroup$
23
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • 6
    $\begingroup$ The RLRsim package is a really convenient way to test random effects using simulation-based likelihood ratio tests. $\endgroup$ – atrichornis Apr 15 '13 at 14:46
  • $\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 '13 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$ – atrichornis Apr 15 '13 at 18:30
  • 5
    $\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 '13 at 13:34
  • 2
    $\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 '13 at 15:16
3
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\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 '13 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 '13 at 8:27
  • $\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 '13 at 12:17
  • 3
    $\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 '13 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$ – Robin Beaumont Jan 21 '18 at 20:10
1
$\begingroup$

I'd like to answer the more 'initial' question.

If you suspect any sort of heterogeneity in variance among either the dependent variable due some some factors, you should go ahead and plot the data using scatter and box plots. Some common patterns to check for, I put this list below from various sources on the web.

Heteroskedasticity Patterns

In addition, plot your dependent variable by factor/treatment groups to see if there is constant variance. If not, you might want to explore random effects or weighted regressions. For eg. this chart below is an example of a funnel shaped variance in my treatment groups. So I choose to go random effects and test the effects on slope and intercepts.

Boxplot to check heteroskedasticity

From here, the answers above cater to your main question. There are also tests that check for heteroskedasticity, one such is here - https://dergipark.org.tr/download/article-file/94971. But i'm not sure if any tests exist to detect group level heteroskedasticity.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Please only use the "Your Answer" field to provide answers to the OP's question. CV is a strict Q&A site, not a discussion forum. The latter, bolded portion of your post is a new question, not an answer to this question. If you have a new question, click the gray ASK QUESTION at the top & ask it there. Since you're new here, you may want to take our tour, which has information for new users. $\endgroup$ – gung - Reinstate Monica♦ Oct 14 '19 at 19:22

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.